module nlte_fomichev

!
! provides calculation of non-LTE heating rates by Fomichev parameterization
!
  use ppgrid,             only: pcols, pver, pverp
  use shr_kind_mod,       only: r8 => shr_kind_r8
  use physconst,          only: r_universal, rearth, avogad, boltz, pi
  use chem_surfvals,      only: chem_surfvals_get
  use cam_abortutils,     only: endrun
  use phys_grid,          only: get_rlon_p, get_rlat_p
  use cam_logfile,        only: iulog
  use spmd_utils,         only: masterproc

  implicit none
  private
  save

! Public interfaces
  public &
     nlte_fomichev_init, &
     nlte_fomichev_calc, &
     nocooling,          &
     o3pcooling

! Private module data

!
! Fomichev radiation parameters
!
   integer :: nrfmc,nrfmg,nrfm,nrfmnlte,nrfmlte,nrfmlteo3,nrfmltelv,nrfmco2

   parameter (nrfmc=59)           ! no. of levels of Fomichev parameterization (2-16.5 by 0.25)
   parameter (nrfmg=8)            ! no. of levels between ground and first calculated level (0-1.75 by 0.25)
   parameter (nrfm=nrfmc+nrfmg)   ! total no. of levels of Fomichev paramterization
   parameter (nrfmnlte=17)        ! no. of levels of NLTE calculation + 1(b.c. at 12.5)
   parameter (nrfmlte=43)         ! no. of levels of LTE calculation
   parameter (nrfmlteo3=35)       ! no. of levels of LTE calculation - O3 ONLY!
   parameter (nrfmltelv=9)        ! no. of levels used in the LTE integral
   parameter (nrfmco2=4)          ! no. of CO2 precalculated profiles

  integer i,ix,js

  real(r8) :: o1_mw                      ! O molecular weight
  real(r8) :: o2_mw                      ! O2 molecular weight
  real(r8) :: o3_mw                      ! O3 molecular weight
  real(r8) :: co2_mw                     ! CO2 molecular weight
  real(r8) :: n2_mw                      ! N2 molecular weight
  real(r8) :: no_mw                      ! NO molecular weight
  
! Physical constants is cgs units
  real(r8), parameter :: akbl=boltz*1e7_r8   ! Boltzman constant
  real(r8), parameter :: anav=avogad*1e-3_r8 ! Avogadro Number
  real(r8), parameter :: grav0=980._r8       ! gravitational constant
  real(r8)            :: arad                ! planet's radius (cm)
  real(r8), parameter :: ur=r_universal*1e4_r8 ! universal gas constant (R_star)
  
  real(r8), parameter :: a10=1.5988_r8         ! reaction constant
  real(r8), parameter :: const=2.63187E11_r8   ! reaction constant
  real(r8), parameter :: constb=9.08795e9_r8   ! reaction constant

  real(r8), parameter :: ptop_co2cool=7.42e-3_r8 ! top pressure level for co2 cool calculation (Pa)
  integer :: ktop_co2cool                        ! the level index defining the top of CO2 cool calculation

!VERTICAL GRIDs to be used in IR scheme
!XR(67) - pressure scale heights, psh's, (=0-16.5) at which input parameter
!         should be given
  real(r8) xr(nrfm)

!DATA for "LTE" parameterization (x=2-12.5)
!IG(9) - indexes of level which should be used to account for the internal
!        atmospheric heat exchange
!AO3(35,9) - coefficients for O3 scheme to calculate cooling rate in
!            "erg/g/s" at levels x=2-10.5 (with a step of 0.25)
!CO2O(4) - vmr for basic CO2
!"LTE-coefficients" for CO2 scheme using to calculate cooling rate in
!"erg/g/s" in region x=2-12.5 (with a step of 0.25). To account for internal
!heat exchange 9 level in atmosphere are needed. 
! A150, B150(43,9) - for 150ppm of CO2
! A360, B360(43,9) - for 360ppm of CO2
! A540, B540(43,9) - for 540ppm of CO2
! A720, B720(43,9) - for 720ppm of CO2
   integer ig(nrfmltelv)

   real(r8) a150(nrfmlte,nrfmltelv)
   real(r8) b150(nrfmlte,nrfmltelv)
   real(r8) a360(nrfmlte,nrfmltelv)
   real(r8) b360(nrfmlte,nrfmltelv)
   real(r8) a540(nrfmlte,nrfmltelv)
   real(r8) b540(nrfmlte,nrfmltelv)
   real(r8) a720(nrfmlte,nrfmltelv)
   real(r8) b720(nrfmlte,nrfmltelv)
   real(r8) co2o(nrfmco2)

!DATA for NLTE parameterization CO2 (x=12.5-16.5)
!UCO2RO, ALO(51) - CO2 column amount and corresponding escape functions
!                                             (eventually, their log)
!COR150, COR360, COR540, COR720(6) - correction to escape functions to
!  calculate coefficients for the reccurence formula between x=12.5 and 13.75
!UCO2CO(6) - CO2 column amount at x=12.5-13.75 (step - 0.25) for 360 ppm
!            to be used to correct escape functions in this region
   real(r8) uco2ro(51)
   real(r8) alo(51)
   real(r8) cor150(6)
   real(r8) cor360(6)
   real(r8) cor540(6)
   real(r8) cor720(6)
   real(r8) uco2co(6)
   real(r8) ao3(nrfmlteo3,nrfmltelv)

   data (xr(i), i=1,67)/                                                 &
       0.0_r8, 0.25_r8, 0.5_r8, 0.75_r8, 1.0_r8, 1.25_r8, 1.5_r8, 1.75_r8, 2.0_r8, 2.25_r8, 2.5_r8, 2.75_r8, &
       3.0_r8, 3.25_r8, 3.5_r8, 3.75_r8, 4.0_r8, 4.25_r8, 4.5_r8, 4.75_r8, 5.0_r8, 5.25_r8, 5.5_r8, 5.75_r8, &
       6.0_r8, 6.25_r8, 6.5_r8, 6.75_r8, 7.0_r8, 7.25_r8, 7.5_r8, 7.75_r8, 8.0_r8, 8.25_r8, 8.5_r8, 8.75_r8, &
       9.0_r8, 9.25_r8, 9.5_r8, 9.75_r8,10.0_r8,10.25_r8,10.5_r8,10.75_r8,11.0_r8,11.25_r8,11.5_r8,11.75_r8, &
      12.0_r8,12.25_r8,12.5_r8,12.75_r8,13.0_r8,13.25_r8,13.5_r8,13.75_r8,14.0_r8,14.25_r8,14.5_r8,14.75_r8, &
      15.0_r8,15.25_r8,15.5_r8,15.75_r8,16.0_r8,16.25_r8,16.5_r8/

   data (ig(i),i=1,9)/-25,-12,-7,-3,-1,0,1,3,6/


   data (co2o(i),i=1,4)/150.e-6_r8, 360.e-6_r8, 540.e-6_r8, 720.e-6_r8/

   data (uco2co(i),i=1,6) /2.546953e+16_r8,1.913609e+16_r8,1.425730e+16_r8,      &
                              1.052205e+16_r8,7.700499e+15_r8,5.596946e+15_r8/
   data (cor150(i),i=1,6) /2.530798e-01_r8,2.048590e-01_r8,1.050616e-01_r8,      &
                              9.172653e-02_r8,4.528593e-02_r8,3.384089e-02_r8/
   data (cor360(i),i=1,6) /4.848122e-01_r8,4.224252e-01_r8,2.958043e-01_r8,      &
                              2.067767e-01_r8,1.244037e-01_r8,5.792163e-02_r8/
   data (cor540(i),i=1,6) /6.263170e-01_r8,5.332695e-01_r8,3.773434e-01_r8,      &
                              2.592216e-01_r8,1.514242e-01_r8,6.889337e-02_r8/
   data (cor720(i),i=1,6) /7.248522e-01_r8,6.199704e-01_r8,4.525580e-01_r8,      &
                              3.046277e-01_r8,1.787235e-01_r8,7.953181e-02_r8/

   data (uco2ro(i),i=1,51)   /2.699726e+11_r8,5.810773e+11_r8,1.106722e+12_r8,   &
       1.952319e+12_r8,3.306797e+12_r8,5.480155e+12_r8,8.858565e+12_r8,1.390142e+13_r8,&
       2.129301e+13_r8,3.209300e+13_r8,4.784654e+13_r8,7.091442e+13_r8,1.052353e+14_r8,&
       1.565317e+14_r8,2.320320e+14_r8,3.415852e+14_r8,4.986668e+14_r8,7.212717e+14_r8,&
       1.033831e+15_r8,1.469497e+15_r8,2.073209e+15_r8,2.905406e+15_r8,4.044901e+15_r8,&
       5.596946e+15_r8,7.700499e+15_r8,1.052205e+16_r8,1.425730e+16_r8,1.913609e+16_r8,&
       2.546953e+16_r8,3.366464e+16_r8,4.421144e+16_r8,5.775381e+16_r8,7.514254e+16_r8,&
       9.747013e+16_r8,1.261393e+17_r8,1.629513e+17_r8,2.102188e+17_r8,2.709114e+17_r8,&
       3.488423e+17_r8,4.489076e+17_r8,5.773939e+17_r8,7.423736e+17_r8,9.542118e+17_r8,&
       1.226217e+18_r8,1.575480e+18_r8,2.023941e+18_r8,2.599777e+18_r8,3.339164e+18_r8,&
       4.288557e+18_r8,5.507602e+18_r8,7.072886e+18_r8/

   data (alo(i),i=1,51) /-2.410106e-04_r8,-5.471415e-04_r8,-1.061586e-03_r8,     &
              -1.879789e-03_r8,-3.166020e-03_r8,-5.185436e-03_r8,-8.216667e-03_r8,  &
              -1.250894e-02_r8,-1.838597e-02_r8,-2.631114e-02_r8,-3.688185e-02_r8,  &  
              -5.096491e-02_r8,-7.004056e-02_r8,-9.603746e-02_r8,-1.307683e-01_r8,  &
              -1.762946e-01_r8,-2.350226e-01_r8,-3.095215e-01_r8,-4.027339e-01_r8,  &
              -5.178570e-01_r8,-6.581256e-01_r8,-8.265003e-01_r8,-1.024684e+00_r8,  &
              -1.252904e+00_r8,-1.509470e+00_r8,-1.788571e+00_r8,-2.081700e+00_r8,  &
              -2.379480e+00_r8,-2.675720e+00_r8,-2.967325e+00_r8,-3.252122e+00_r8,  &
              -3.530485e+00_r8,-3.803720e+00_r8,-4.072755e+00_r8,-4.338308e+00_r8,  &
              -4.601048e+00_r8,-4.861585e+00_r8,-5.120370e+00_r8,-5.377789e+00_r8,  &
              -5.634115e+00_r8,-5.889388e+00_r8,-6.143488e+00_r8,-6.396436e+00_r8,  &
              -6.648774e+00_r8,-6.901465e+00_r8,-7.155207e+00_r8,-7.409651e+00_r8,  &
              -7.663536e+00_r8,-7.915682e+00_r8,-8.165871e+00_r8,-8.415016e+00_r8/

   data ((ao3(ix,js),js=1,9),ix=1,10)/                                  &
       5.690e+09_r8, 0.000e+00_r8, 1.962e+09_r8, 5.015e+09_r8, 7.105e+09_r8,-3.952e+10_r8,&
       8.916e+09_r8, 1.204e+09_r8, 5.434e+09_r8, 4.997e+09_r8, 0.000e+00_r8, 2.168e+09_r8,&
       4.981e+09_r8, 7.327e+09_r8,-3.855e+10_r8, 9.296e+09_r8, 1.326e+09_r8, 4.575e+09_r8,&
       4.179e+09_r8, 0.000e+00_r8, 2.198e+09_r8, 4.594e+09_r8, 7.780e+09_r8,-3.689e+10_r8,&
       9.624e+09_r8, 1.107e+09_r8, 3.736e+09_r8, 3.469e+09_r8, 0.000e+00_r8, 1.941e+09_r8,&
       3.921e+09_r8, 8.423e+09_r8,-3.508e+10_r8, 9.859e+09_r8, 8.736e+08_r8, 3.058e+09_r8,&
       2.930e+09_r8, 0.000e+00_r8, 1.649e+09_r8, 3.378e+09_r8, 8.985e+09_r8,-3.362e+10_r8,&
       9.941e+09_r8, 8.110e+08_r8, 2.367e+09_r8, 2.020e+09_r8, 3.866e+08_r8, 1.368e+09_r8,&
       3.049e+09_r8, 9.452e+09_r8,-3.249e+10_r8, 1.010e+10_r8, 6.212e+08_r8, 2.117e+09_r8,&
       1.652e+09_r8, 3.711e+08_r8, 1.175e+09_r8, 2.747e+09_r8, 9.818e+09_r8,-3.165e+10_r8,&
       1.014e+10_r8, 5.470e+08_r8, 1.843e+09_r8, 1.391e+09_r8, 3.497e+08_r8, 1.033e+09_r8,&
       2.570e+09_r8, 1.008e+10_r8,-3.102e+10_r8, 1.017e+10_r8, 5.121e+08_r8, 1.552e+09_r8,&
       1.197e+09_r8, 3.329e+08_r8, 9.689e+08_r8, 2.427e+09_r8, 1.024e+10_r8,-3.056e+10_r8,&
       1.011e+10_r8, 4.521e+08_r8, 1.312e+09_r8, 9.961e+08_r8, 3.390e+08_r8, 9.867e+08_r8,&
       2.375e+09_r8, 1.032e+10_r8,-3.033e+10_r8, 1.011e+10_r8, 5.322e+08_r8, 1.178e+09_r8/

   data ((ao3(ix,js),js=1,9),ix=11,20)/                                 &
       8.408e+08_r8, 4.922e+08_r8, 9.226e+08_r8, 2.358e+09_r8, 1.039e+10_r8,-3.030e+10_r8,&
       1.004e+10_r8, 5.551e+08_r8, 1.141e+09_r8, 7.250e+08_r8, 5.353e+08_r8, 9.493e+08_r8,&
       2.497e+09_r8, 1.040e+10_r8,-3.047e+10_r8, 9.971e+09_r8, 8.389e+08_r8, 9.122e+08_r8,&
       6.484e+08_r8, 5.214e+08_r8, 1.014e+09_r8, 2.625e+09_r8, 1.038e+10_r8,-3.058e+10_r8,&
       9.903e+09_r8, 7.720e+08_r8, 9.455e+08_r8, 5.785e+08_r8, 4.987e+08_r8, 1.037e+09_r8,&
       2.694e+09_r8, 1.039e+10_r8,-3.089e+10_r8, 9.698e+09_r8, 1.068e+09_r8, 8.903e+08_r8,&
       5.272e+08_r8, 4.955e+08_r8, 1.194e+09_r8, 2.998e+09_r8, 1.026e+10_r8,-3.160e+10_r8,&
       9.466e+09_r8, 1.239e+09_r8, 9.659e+08_r8, 4.891e+08_r8, 5.131e+08_r8, 1.299e+09_r8,&
       3.456e+09_r8, 1.016e+10_r8,-3.275e+10_r8, 9.003e+09_r8, 1.593e+09_r8, 1.083e+09_r8,&
       4.643e+08_r8, 5.668e+08_r8, 1.500e+09_r8, 4.096e+09_r8, 9.972e+09_r8,-3.422e+10_r8,&
       8.345e+09_r8, 1.790e+09_r8, 1.256e+09_r8, 4.461e+08_r8, 6.391e+08_r8, 1.870e+09_r8,&
       4.881e+09_r8, 9.608e+09_r8,-3.621e+10_r8, 7.386e+09_r8, 2.105e+09_r8, 1.460e+09_r8,&
       6.471e+08_r8, 1.536e+08_r8, 2.567e+09_r8, 5.909e+09_r8, 9.088e+09_r8,-3.844e+10_r8,&
       6.248e+09_r8, 2.173e+09_r8, 1.536e+09_r8, 5.377e+08_r8, 1.101e+08_r8, 3.474e+09_r8,&
       6.992e+09_r8, 8.333e+09_r8,-4.065e+10_r8, 4.979e+09_r8, 2.019e+09_r8, 1.594e+09_r8/

   data ((ao3(ix,js),js=1,9),ix=21,30)/                                 &
       7.064e+08_r8, 1.346e+08_r8, 4.535e+09_r8, 7.814e+09_r8, 7.361e+09_r8,-4.248e+10_r8,&
       3.786e+09_r8, 1.733e+09_r8, 1.431e+09_r8, 1.210e+09_r8, 2.133e+08_r8, 5.772e+09_r8,&
       8.402e+09_r8, 6.184e+09_r8,-4.400e+10_r8, 2.883e+09_r8, 1.530e+09_r8, 1.163e+09_r8,&
       2.121e+09_r8, 3.450e+08_r8, 7.111e+09_r8, 8.611e+09_r8, 4.874e+09_r8,-4.499e+10_r8,&
       2.055e+09_r8, 1.126e+09_r8, 8.840e+08_r8, 3.393e+09_r8, 5.066e+08_r8, 8.176e+09_r8,&
       8.419e+09_r8, 3.777e+09_r8,-4.573e+10_r8, 1.418e+09_r8, 8.927e+08_r8, 6.356e+08_r8,&
       4.838e+09_r8, 6.534e+08_r8, 9.188e+09_r8, 7.895e+09_r8, 2.746e+09_r8,-4.628e+10_r8,&
       9.548e+08_r8, 6.033e+08_r8, 4.810e+08_r8, 5.852e+09_r8, 1.286e+09_r8, 9.867e+09_r8,&
       7.121e+09_r8, 2.046e+09_r8,-4.663e+10_r8, 6.643e+08_r8, 4.034e+08_r8, 2.258e+08_r8,&
       6.624e+09_r8, 2.439e+09_r8, 1.044e+10_r8, 5.950e+09_r8, 1.380e+09_r8,-4.691e+10_r8,&
       4.209e+08_r8, 2.635e+08_r8, 2.211e+08_r8, 7.030e+09_r8, 3.865e+09_r8, 1.055e+10_r8,&
       4.781e+09_r8, 1.009e+09_r8,-4.708e+10_r8, 3.238e+08_r8, 2.107e+08_r8, 1.094e+08_r8,&
       8.099e+09_r8, 4.843e+09_r8, 1.065e+10_r8, 3.732e+09_r8, 6.798e+08_r8,-4.719e+10_r8,&
       1.852e+08_r8, 1.233e+08_r8, 1.513e+08_r8, 8.382e+09_r8, 6.426e+09_r8, 1.018e+10_r8,&
       2.771e+09_r8, 4.593e+08_r8,-4.726e+10_r8, 1.084e+08_r8, 8.368e+07_r8, 0.000e+00_r8/

   data ((ao3(ix,js),js=1,9),ix=31,35)/                                 &
       9.545e+09_r8, 7.667e+09_r8, 9.429e+09_r8, 1.927e+09_r8, 3.143e+08_r8,-4.731e+10_r8,&
       7.075e+07_r8, 5.179e+07_r8, 0.000e+00_r8, 1.036e+10_r8, 9.146e+09_r8, 8.014e+09_r8,&
       1.355e+09_r8, 2.280e+08_r8,-4.734e+10_r8, 4.459e+07_r8, 3.463e+07_r8, 0.000e+00_r8,&
       1.180e+10_r8, 1.021e+10_r8, 6.725e+09_r8, 9.575e+08_r8, 1.091e+08_r8,-4.736e+10_r8,&
       3.018e+07_r8, 2.000e+07_r8, 0.000e+00_r8, 8.878e+09_r8, 7.184e+09_r8, 3.569e+09_r8,&
       4.023e+08_r8, 4.656e+07_r8,-3.079e+10_r8, 1.374e+07_r8, 0.000e+00_r8, 0.000e+00_r8,&
       4.652e+09_r8, 3.469e+09_r8, 1.385e+09_r8, 1.114e+08_r8, 1.368e+07_r8,-1.421e+10_r8,&
       3.553e+06_r8, 0.000e+00_r8, 0.000e+00_r8/

   data ((a150(ix,js),js=1,9),ix=1,5)/                                  &      
       2.5035e+01_r8, 0.0000e+00_r8, 9.4137e+01_r8, 2.4765e+03_r8, 3.2467e+03_r8,      &
      -1.1090e+04_r8, 4.0701e+02_r8, 2.1062e+03_r8, 6.1997e+02_r8, 3.1497e+01_r8,      &
       0.0000e+00_r8, 1.5661e+02_r8, 3.1886e+03_r8, 3.0556e+03_r8,-1.2560e+04_r8,      &
       2.1713e+02_r8, 2.2157e+03_r8, 7.7184e+02_r8, 3.6209e+01_r8, 0.0000e+00_r8,      &
       2.5391e+02_r8, 3.8867e+03_r8, 2.8525e+03_r8,-1.4009e+04_r8,-3.8903e+00_r8,      &
       2.2576e+03_r8, 9.1027e+02_r8, 3.9456e+01_r8, 0.0000e+00_r8, 4.0539e+02_r8,      &
       4.4674e+03_r8, 2.6199e+03_r8,-1.5203e+04_r8,-2.0232e+02_r8, 2.2454e+03_r8,      &
       1.0123e+03_r8, 4.3855e+01_r8, 0.0000e+00_r8, 6.3518e+02_r8, 4.9490e+03_r8,      &
       2.3873e+03_r8,-1.6282e+04_r8,-4.0148e+02_r8, 2.2334e+03_r8, 1.0983e+03_r8/ 

   data ((a150(ix,js),js=1,9),ix=6,10)/                                 &
       1.4880e+01_r8, 5.2053e+01_r8, 8.8187e+02_r8, 5.4055e+03_r8, 2.1311e+03_r8,      &
      -1.7305e+04_r8,-5.7150e+02_r8, 2.2465e+03_r8, 1.1736e+03_r8, 1.8779e+01_r8,      &
       8.2377e+01_r8, 1.1180e+03_r8, 5.7943e+03_r8, 1.8908e+03_r8,-1.8247e+04_r8,      &
      -7.1038e+02_r8, 2.3030e+03_r8, 1.2335e+03_r8, 2.3004e+01_r8, 1.2959e+02_r8,      &
       1.3044e+03_r8, 6.1369e+03_r8, 1.6879e+03_r8,-1.9121e+04_r8,-8.4776e+02_r8,      &
       2.4019e+03_r8, 1.2764e+03_r8, 2.7051e+01_r8, 1.9749e+02_r8, 1.4427e+03_r8,      &
       6.4390e+03_r8, 1.4997e+03_r8,-1.9896e+04_r8,-9.8529e+02_r8, 2.5405e+03_r8,      &
       1.2946e+03_r8, 3.0107e+01_r8, 2.8477e+02_r8, 1.5605e+03_r8, 6.6846e+03_r8,      &
       1.3469e+03_r8,-2.0582e+04_r8,-1.1115e+03_r8, 2.7010e+03_r8, 1.2951e+03_r8/

   data ((a150(ix,js),js=1,9),ix=11,15)/                                &
       3.3100e+01_r8, 3.8721e+02_r8, 1.6772e+03_r8, 6.8908e+03_r8, 1.2666e+03_r8,      &
      -2.1299e+04_r8,-1.2333e+03_r8, 2.9002e+03_r8, 1.2922e+03_r8, 3.6769e+01_r8,      &
       5.0026e+02_r8, 1.8186e+03_r8, 7.0844e+03_r8, 1.2377e+03_r8,-2.2143e+04_r8,      &
      -1.3997e+03_r8, 3.1384e+03_r8, 1.3038e+03_r8, 4.0958e+01_r8, 6.2580e+02_r8,      &
       1.9791e+03_r8, 7.3078e+03_r8, 1.1888e+03_r8,-2.3138e+04_r8,-1.5839e+03_r8,      &
       3.4194e+03_r8, 1.3343e+03_r8, 4.5506e+01_r8, 7.5175e+02_r8, 2.1209e+03_r8,      &
       7.5851e+03_r8, 1.1225e+03_r8,-2.4221e+04_r8,-1.8076e+03_r8, 3.7350e+03_r8,      &
       1.3778e+03_r8, 5.0606e+01_r8, 8.8360e+02_r8, 2.2346e+03_r8, 7.9275e+03_r8,      &
       1.0078e+03_r8,-2.5290e+04_r8,-2.0653e+03_r8, 4.0577e+03_r8, 1.4264e+03_r8/

   data ((a150(ix,js),js=1,9),ix=16,20)/                                &
       5.6331e+01_r8, 1.0256e+03_r8, 2.3251e+03_r8, 8.3240e+03_r8, 8.5705e+02_r8,      &
      -2.6376e+04_r8,-2.3440e+03_r8, 4.3896e+03_r8, 1.4827e+03_r8, 6.2586e+01_r8,      &
       1.1799e+03_r8, 2.4192e+03_r8, 8.7537e+03_r8, 7.1111e+02_r8,-2.7532e+04_r8,      &
      -2.6284e+03_r8, 4.7408e+03_r8, 1.5489e+03_r8, 6.9499e+01_r8, 1.3513e+03_r8,      &
       2.5385e+03_r8, 9.2024e+03_r8, 5.9855e+02_r8,-2.8837e+04_r8,-2.9066e+03_r8,      &
       5.1230e+03_r8, 1.6315e+03_r8, 9.8718e+01_r8, 1.4643e+03_r8, 2.7198e+03_r8,      &
       9.6727e+03_r8, 5.6542e+02_r8,-3.0324e+04_r8,-3.1708e+03_r8, 5.5386e+03_r8,      &
       1.7411e+03_r8, 1.4214e+02_r8, 1.5587e+03_r8, 2.9412e+03_r8, 1.0183e+04_r8,      &
       6.2959e+02_r8,-3.2074e+04_r8,-3.4238e+03_r8, 5.9832e+03_r8, 1.8906e+03_r8/

   data ((a150(ix,js),js=1,9),ix=21,25)/                                &
       2.0882e+02_r8, 1.6315e+03_r8, 3.1941e+03_r8, 1.0757e+04_r8, 8.0531e+02_r8,      &
      -3.4141e+04_r8,-3.6662e+03_r8, 6.4368e+03_r8, 2.0964e+03_r8, 3.0513e+02_r8,      &
       1.6987e+03_r8, 3.4635e+03_r8, 1.1429e+04_r8, 1.0747e+03_r8,-3.6585e+04_r8,      &
      -3.8843e+03_r8, 6.8932e+03_r8, 2.3807e+03_r8, 4.3086e+02_r8, 1.7634e+03_r8,      &
       3.7650e+03_r8, 1.2253e+04_r8, 1.4449e+03_r8,-3.9487e+04_r8,-4.0949e+03_r8,      &
       7.3521e+03_r8, 2.7604e+03_r8, 5.6228e+02_r8, 1.8545e+03_r8, 4.0922e+03_r8,      &
       1.3286e+04_r8, 1.8655e+03_r8,-4.2969e+04_r8,-4.3141e+03_r8, 7.7726e+03_r8,      &
       3.3137e+03_r8, 6.5591e+02_r8, 1.9789e+03_r8, 4.4192e+03_r8, 1.4473e+04_r8,      &
       2.2732e+03_r8,-4.6758e+04_r8,-4.5119e+03_r8, 8.1060e+03_r8, 4.0069e+03_r8/

   data ((a150(ix,js),js=1,9),ix=26,31)/                                &
       7.1489e+02_r8, 2.1391e+03_r8, 4.7555e+03_r8, 1.5792e+04_r8, 2.6329e+03_r8,      &
      -5.0784e+04_r8,-4.6936e+03_r8, 8.3781e+03_r8, 4.7906e+03_r8, 7.5917e+02_r8,      &
       2.3306e+03_r8, 5.1321e+03_r8, 1.7236e+04_r8, 2.9779e+03_r8,-5.5179e+04_r8,      &
      -4.8828e+03_r8, 8.7136e+03_r8, 5.5138e+03_r8, 8.1061e+02_r8, 2.5462e+03_r8,      &
       5.5725e+03_r8, 1.8761e+04_r8, 3.3368e+03_r8,-6.0007e+04_r8,-5.0788e+03_r8,      &
       9.1950e+03_r8, 6.1026e+03_r8, 8.6939e+02_r8, 2.7684e+03_r8, 6.1037e+03_r8,      &
       2.0165e+04_r8, 3.7711e+03_r8,-6.4792e+04_r8,-5.2866e+03_r8, 9.9049e+03_r8,      &
       6.1727e+03_r8, 9.4294e+02_r8, 3.0123e+03_r8, 6.6875e+03_r8, 2.1225e+04_r8,      &
       4.4266e+03_r8,-6.9341e+04_r8,-5.5392e+03_r8, 1.1040e+04_r8, 5.5220e+03_r8,      &
       1.0189e+03_r8, 3.2497e+03_r8, 7.2300e+03_r8, 2.0698e+04_r8, 5.1305e+03_r8,      &
      -7.1191e+04_r8,-5.7478e+03_r8, 1.2510e+04_r8, 3.5534e+03_r8/

   data ((a150(ix,js),js=1,9),ix=32,37)/                                &
       1.1449e+03_r8, 3.7084e+03_r8, 7.4630e+03_r8, 1.9495e+04_r8, 5.5511e+03_r8,      &
      -7.3819e+04_r8,-5.4997e+03_r8, 1.4972e+04_r8, 2.3195e+03_r8, 1.2999e+03_r8,      &
       4.3325e+03_r8, 7.4755e+03_r8, 1.6338e+04_r8, 5.1750e+03_r8,-7.5166e+04_r8,      &
      -4.4753e+03_r8, 1.7530e+04_r8, 2.4555e+03_r8, 1.4512e+03_r8, 5.0331e+03_r8,      &
       7.0224e+03_r8, 1.5180e+04_r8, 3.2675e+03_r8,-7.8134e+04_r8,-3.0385e+03_r8,      &
       2.0430e+04_r8, 2.7439e+03_r8, 1.5838e+03_r8, 5.7655e+03_r8, 6.0092e+03_r8,      &
       1.4199e+04_r8, 5.0935e+03_r8,-8.0938e+04_r8,-3.1873e+03_r8, 2.1397e+04_r8,      &
       2.3826e+03_r8, 1.6805e+03_r8, 6.6768e+03_r8, 4.2956e+03_r8, 1.4811e+04_r8,      &
       8.4327e+03_r8,-8.6256e+04_r8,-4.1912e+03_r8, 2.2919e+04_r8, 1.4885e+03_r8,      &
       1.7701e+03_r8, 7.7803e+03_r8, 3.4348e+03_r8, 1.8569e+04_r8, 1.3687e+04_r8,      &
      -1.0164e+05_r8,-6.1459e+03_r8, 2.7017e+04_r8, 1.4367e+03_r8/

   data ((a150(ix,js),js=1,9),ix=38,43)/                                &
       2.1081e+03_r8, 9.8559e+03_r8, 3.3385e+03_r8, 2.5453e+04_r8, 1.9163e+04_r8,      &
      -1.3508e+05_r8,-7.9100e+03_r8, 4.0371e+04_r8, 3.1637e+03_r8, 2.5117e+03_r8,      &
       1.2043e+04_r8, 2.9837e+03_r8, 3.0525e+04_r8, 2.1190e+04_r8,-1.6547e+05_r8,      &
      -1.0361e+04_r8, 5.5753e+04_r8, 6.3850e+03_r8, 2.7271e+03_r8, 1.3131e+04_r8,      &
       1.4005e+03_r8, 3.1777e+04_r8, 1.8261e+04_r8,-1.7548e+05_r8,-1.2891e+04_r8,      &
       6.4744e+04_r8, 9.9029e+03_r8, 2.4817e+03_r8, 1.2549e+04_r8, 1.2521e+03_r8,      &
       3.2506e+04_r8, 1.5901e+04_r8,-1.6344e+05_r8,-1.7080e+04_r8, 5.9506e+04_r8,      &
       7.9804e+03_r8, 2.0450e+03_r8, 1.1403e+04_r8, 3.8851e+03_r8, 3.7098e+04_r8,      &
       1.3215e+04_r8,-1.5453e+05_r8,-2.1257e+04_r8, 5.2256e+04_r8, 5.4705e+03_r8,      &
       1.7609e+03_r8, 1.1782e+04_r8, 7.3394e+03_r8, 3.6003e+04_r8, 5.7292e+03_r8,      &
      -1.4489e+05_r8,-2.1518e+04_r8, 4.7419e+04_r8, 2.0765e+03_r8/

   data ((b150(ix,js),js=1,9),ix=1,5)/                                  &
       3.5399e+03_r8, 0.0000e+00_r8, 6.0508e+03_r8, 5.9109e+04_r8, 2.4118e+04_r8,      &
      -1.8103e+05_r8,-7.3845e+02_r8, 1.8008e+04_r8, 1.7036e+04_r8, 3.1195e+03_r8,      &
       0.0000e+00_r8, 6.9773e+03_r8, 5.7201e+04_r8, 2.0239e+04_r8,-1.7100e+05_r8,      &
      -2.3973e+03_r8, 1.8321e+04_r8, 1.5737e+04_r8, 3.0408e+03_r8, 0.0000e+00_r8,      &
       9.0810e+03_r8, 5.5104e+04_r8, 1.3930e+04_r8,-1.6106e+05_r8,-4.9406e+03_r8,      &
       1.6082e+04_r8, 1.4928e+04_r8, 3.2990e+03_r8, 0.0000e+00_r8, 1.3161e+04_r8,      &
       5.6155e+04_r8, 1.0436e+04_r8,-1.6715e+05_r8,-7.0673e+03_r8, 1.5202e+04_r8,      &
       1.5529e+04_r8, 3.5043e+03_r8, 0.0000e+00_r8, 1.8035e+04_r8, 5.4916e+04_r8,      &
       7.8353e+03_r8,-1.7025e+05_r8,-8.4919e+03_r8, 1.4317e+04_r8, 1.5764e+04_r8/

   data ((b150(ix,js),js=1,9),ix=6,10)/                                 &
       1.8653e+03_r8, 2.7059e+03_r8, 1.9745e+04_r8, 5.4113e+04_r8, 5.9014e+03_r8,      &
      -1.7121e+05_r8,-9.3646e+03_r8, 1.3888e+04_r8, 1.5667e+04_r8, 2.0888e+03_r8,      &
       3.6216e+03_r8, 1.9181e+04_r8, 5.2430e+04_r8, 4.3250e+03_r8,-1.6730e+05_r8,      &
      -9.8248e+03_r8, 1.3662e+04_r8, 1.5012e+04_r8, 2.2582e+03_r8, 4.7714e+03_r8,      &
       1.7132e+04_r8, 4.9716e+04_r8, 2.5383e+03_r8,-1.5818e+05_r8,-1.0089e+04_r8,      &
       1.3299e+04_r8, 1.3889e+04_r8, 2.4496e+03_r8, 6.2997e+03_r8, 1.5321e+04_r8,      &
       4.8636e+04_r8, 9.3940e+02_r8,-1.5318e+05_r8,-1.0680e+04_r8, 1.3768e+04_r8,      &
       1.2975e+04_r8, 2.7048e+03_r8, 8.2966e+03_r8, 1.4669e+04_r8, 4.9707e+04_r8,      &
      -6.8884e+02_r8,-1.5464e+05_r8,-1.1834e+04_r8, 1.4872e+04_r8, 1.2636e+04_r8/

   data ((b150(ix,js),js=1,9),ix=11,15)/                                &
       3.0080e+03_r8, 1.0337e+04_r8, 1.4758e+04_r8, 5.1024e+04_r8,-2.2171e+03_r8,      &
      -1.5788e+05_r8,-1.3463e+04_r8, 1.6569e+04_r8, 1.2308e+04_r8, 3.3257e+03_r8,      &
       1.2025e+04_r8, 1.5122e+04_r8, 5.1259e+04_r8,-3.7550e+03_r8,-1.5938e+05_r8,      &
      -1.5169e+04_r8, 1.8104e+04_r8, 1.1931e+04_r8, 3.6784e+03_r8, 1.3370e+04_r8,      &
       1.5491e+04_r8, 5.0380e+04_r8,-5.6000e+03_r8,-1.5854e+05_r8,-1.6538e+04_r8,      &
       1.9311e+04_r8, 1.1575e+04_r8, 4.0171e+03_r8, 1.4413e+04_r8, 1.5479e+04_r8,      &
       4.8944e+04_r8,-7.1724e+03_r8,-1.5614e+05_r8,-1.7640e+04_r8, 2.0517e+04_r8,      &
       1.1214e+04_r8, 4.3161e+03_r8, 1.5290e+04_r8, 1.5284e+04_r8, 4.8663e+04_r8,      &
      -7.7645e+03_r8,-1.5691e+05_r8,-1.8632e+04_r8, 2.2492e+04_r8, 1.1169e+04_r8/

   data ((b150(ix,js),js=1,9),ix=16,20)/                                &
       4.6114e+03_r8, 1.6154e+04_r8, 1.4979e+04_r8, 4.9364e+04_r8,-7.9418e+03_r8,      &
      -1.6022e+05_r8,-1.9786e+04_r8, 2.5014e+04_r8, 1.1385e+04_r8, 4.8886e+03_r8,      &
       1.7121e+04_r8, 1.4681e+04_r8, 5.0603e+04_r8,-8.0295e+03_r8,-1.6487e+05_r8,      &
      -2.1104e+04_r8, 2.7565e+04_r8, 1.1912e+04_r8, 5.1220e+03_r8, 1.8151e+04_r8,      &
       1.4554e+04_r8, 5.2649e+04_r8,-8.0370e+03_r8,-1.7162e+05_r8,-2.2730e+04_r8,      &
       3.0243e+04_r8, 1.2705e+04_r8, 6.9450e+03_r8, 1.4807e+04_r8, 1.6724e+04_r8,      &
       5.6069e+04_r8,-8.0590e+03_r8,-1.8189e+05_r8,-2.4936e+04_r8, 3.3115e+04_r8,      &
       1.4017e+04_r8, 9.4904e+03_r8, 1.0887e+04_r8, 1.9534e+04_r8, 6.0737e+04_r8,      &
      -8.1973e+03_r8,-1.9502e+05_r8,-2.7754e+04_r8, 3.5773e+04_r8, 1.5904e+04_r8/

   data ((b150(ix,js),js=1,9),ix=21,25)/                                &
       1.3201e+04_r8, 6.7672e+03_r8, 2.2982e+04_r8, 6.6827e+04_r8,-8.7998e+03_r8,      &
      -2.1116e+05_r8,-3.1231e+04_r8, 3.7848e+04_r8, 1.8559e+04_r8, 1.8353e+04_r8,      &
       2.9582e+03_r8, 2.7018e+04_r8, 7.4757e+04_r8,-1.0238e+04_r8,-2.3087e+05_r8,      &
      -3.5386e+04_r8, 3.9278e+04_r8, 2.2035e+04_r8, 2.4236e+04_r8,-7.1673e+00_r8,      &
       3.1406e+04_r8, 8.4123e+04_r8,-1.2597e+04_r8,-2.5210e+05_r8,-3.9944e+04_r8,      &
       3.9900e+04_r8, 2.5639e+04_r8, 2.8899e+04_r8,-1.8050e+03_r8, 3.5598e+04_r8,      &
       9.3624e+04_r8,-1.6264e+04_r8,-2.7076e+05_r8,-4.3916e+04_r8, 3.8269e+04_r8,      &
       2.9577e+04_r8, 3.2228e+04_r8,-2.3769e+03_r8, 4.1189e+04_r8, 1.0710e+05_r8,      &
      -2.1940e+04_r8,-2.9820e+05_r8,-4.9133e+04_r8, 3.6912e+04_r8, 3.4890e+04_r8/

   data ((b150(ix,js),js=1,9),ix=26,31)/                                &
       3.4901e+04_r8,-1.7054e+03_r8, 4.9616e+04_r8, 1.2729e+05_r8,-2.9857e+04_r8,      &
      -3.4296e+05_r8,-5.6738e+04_r8, 3.7392e+04_r8, 4.2347e+04_r8, 3.7694e+04_r8,      &
       2.6968e+02_r8, 6.1729e+04_r8, 1.5470e+05_r8,-4.0354e+04_r8,-4.0768e+05_r8,      &
      -6.7248e+04_r8, 4.0115e+04_r8, 5.1589e+04_r8, 4.0809e+04_r8, 4.1844e+03_r8,      &
       7.8173e+04_r8, 1.8915e+05_r8,-5.3538e+04_r8,-4.9425e+05_r8,-8.0910e+04_r8,      &
       4.5420e+04_r8, 6.2215e+04_r8, 4.3379e+04_r8, 1.0525e+04_r8, 9.9171e+04_r8,      &
       2.2658e+05_r8,-6.7757e+04_r8,-5.9481e+05_r8,-9.6665e+04_r8, 5.4277e+04_r8,      &
       6.9248e+04_r8, 4.5667e+04_r8, 1.9161e+04_r8, 1.2399e+05_r8, 2.5979e+05_r8,      & 
      -8.0920e+04_r8,-6.9722e+05_r8,-1.1391e+05_r8, 6.8928e+04_r8, 6.7840e+04_r8,      &
       4.8149e+04_r8, 2.9647e+04_r8, 1.5197e+05_r8, 2.6905e+05_r8,-9.0125e+04_r8,      &
      -7.6851e+05_r8,-1.3001e+05_r8, 8.9270e+04_r8, 4.9048e+04_r8/

   data ((b150(ix,js),js=1,9),ix=32,37)/                                &
       5.1966e+04_r8, 4.5702e+04_r8, 1.7720e+05_r8, 2.6058e+05_r8,-9.9247e+04_r8,      &
      -8.4400e+05_r8,-1.4248e+05_r8, 1.1906e+05_r8, 3.5837e+04_r8, 5.7382e+04_r8,      &
       6.7581e+04_r8, 1.9871e+05_r8, 2.1424e+05_r8,-1.0592e+05_r8,-9.0196e+05_r8,      &
      -1.4324e+05_r8, 1.4835e+05_r8, 3.8639e+04_r8, 6.3804e+04_r8, 9.4177e+04_r8,      &
       2.0441e+05_r8, 1.9644e+05_r8,-1.3089e+05_r8,-9.7326e+05_r8,-1.3977e+05_r8,      &
       1.7913e+05_r8, 4.3016e+04_r8, 7.0638e+04_r8, 1.2395e+05_r8, 1.8916e+05_r8,      &
       1.7195e+05_r8,-1.1044e+05_r8,-1.0265e+06_r8,-1.4544e+05_r8, 1.9545e+05_r8,      &
       3.6843e+04_r8, 7.6839e+04_r8, 1.6172e+05_r8, 1.4909e+05_r8, 1.6911e+05_r8,      &
      -7.5704e+04_r8,-1.0973e+06_r8,-1.5936e+05_r8, 2.1662e+05_r8, 2.4064e+04_r8,      &
       8.3232e+04_r8, 2.0819e+05_r8, 1.2835e+05_r8, 2.0677e+05_r8,-3.6307e+04_r8,      &
      -1.2961e+06_r8,-1.9272e+05_r8, 2.6288e+05_r8, 2.1977e+04_r8/
   data ((b150(ix,js),js=1,9),ix=38,43)/                                &
       9.6157e+04_r8, 2.8562e+05_r8, 1.1966e+05_r8, 2.8749e+05_r8,-1.2237e+04_r8,      &
      -1.7528e+06_r8,-2.4263e+05_r8, 4.0970e+05_r8, 4.1317e+04_r8, 1.1357e+05_r8,      &
       3.7596e+05_r8, 9.8652e+04_r8, 3.6598e+05_r8, 1.6096e+04_r8,-2.2962e+06_r8,      &
      -2.8607e+05_r8, 6.2929e+05_r8, 7.4848e+04_r8, 1.2829e+05_r8, 4.4149e+05_r8,      &
       4.2542e+04_r8, 4.3449e+05_r8, 3.5764e+04_r8,-2.7327e+06_r8,-3.1740e+05_r8,      &
       8.5228e+05_r8, 1.2250e+05_r8, 1.2481e+05_r8, 4.3349e+05_r8, 3.3622e+04_r8,      &
       4.8685e+05_r8, 7.3970e+04_r8,-2.7619e+06_r8,-3.5847e+05_r8, 8.7436e+05_r8,      &
       1.0368e+05_r8, 1.0994e+05_r8, 3.7480e+05_r8, 1.2328e+05_r8, 6.1220e+05_r8,      &
       7.4486e+04_r8,-2.8231e+06_r8,-4.3841e+05_r8, 8.4803e+05_r8, 7.9803e+04_r8,      &
       1.0138e+05_r8, 3.7276e+05_r8, 2.4376e+05_r8, 7.1193e+05_r8, 1.4092e+05_r8,      &
      -3.3243e+06_r8,-4.4325e+05_r8, 1.0532e+06_r8, 3.2305e+04_r8/

   data ((a360(ix,js),js=1,9),ix=1,5)/                                  &
       1.7094e+01_r8, 0.0000e+00_r8, 8.2972e+01_r8, 2.5829e+03_r8, 4.4982e+03_r8,      &
      -1.3393e+04_r8, 9.4326e+02_r8, 2.8560e+03_r8, 5.4669e+02_r8, 2.2911e+01_r8,      &
       0.0000e+00_r8, 1.4184e+02_r8, 3.4315e+03_r8, 4.6065e+03_r8,-1.5742e+04_r8,      &
       7.9383e+02_r8, 3.2302e+03_r8, 7.5442e+02_r8, 2.7345e+01_r8, 0.0000e+00_r8,      &
       2.3163e+02_r8, 4.3474e+03_r8, 4.7106e+03_r8,-1.8243e+04_r8, 5.5670e+02_r8,      &
       3.4688e+03_r8, 9.8490e+02_r8, 3.0814e+01_r8, 0.0000e+00_r8, 3.7251e+02_r8,      &
       5.2547e+03_r8, 4.7045e+03_r8,-2.0572e+04_r8, 3.1055e+02_r8, 3.5462e+03_r8,      &
       1.2007e+03_r8, 3.5494e+01_r8, 0.0000e+00_r8, 5.9797e+02_r8, 6.1381e+03_r8,      &
       4.5752e+03_r8,-2.2700e+04_r8,-9.7490e-02_r8, 3.5437e+03_r8, 1.3996e+03_r8/

   data ((a360(ix,js),js=1,9),ix=6,10)/                                 &
       9.1340e+00_r8, 4.8891e+01_r8, 8.6264e+02_r8, 7.0522e+03_r8, 4.2702e+03_r8,      &
      -2.4662e+04_r8,-3.0600e+02_r8, 3.5640e+03_r8, 1.5687e+03_r8, 1.2880e+01_r8,      &
       8.0570e+01_r8, 1.1535e+03_r8, 7.8722e+03_r8, 3.8761e+03_r8,-2.6396e+04_r8,      &
      -5.6181e+02_r8, 3.6560e+03_r8, 1.7011e+03_r8, 1.7276e+01_r8, 1.3152e+02_r8,      &
       1.4156e+03_r8, 8.5681e+03_r8, 3.5175e+03_r8,-2.7925e+04_r8,-8.2205e+02_r8,      &
       3.8262e+03_r8, 1.7943e+03_r8, 2.1698e+01_r8, 2.0535e+02_r8, 1.6357e+03_r8,      &
       9.1928e+03_r8, 3.1758e+03_r8,-2.9308e+04_r8,-1.0872e+03_r8, 4.0624e+03_r8,      &
       1.8441e+03_r8, 2.5396e+01_r8, 2.9819e+02_r8, 1.8384e+03_r8, 9.7178e+03_r8,      &
       2.8742e+03_r8,-3.0523e+04_r8,-1.3368e+03_r8, 4.3248e+03_r8, 1.8525e+03_r8/

   data ((a360(ix,js),js=1,9),ix=11,15)/                                &
       2.8749e+01_r8, 4.0334e+02_r8, 2.0531e+03_r8, 1.0151e+04_r8, 2.6704e+03_r8,      &
      -3.1682e+04_r8,-1.5687e+03_r8, 4.6207e+03_r8, 1.8353e+03_r8, 3.2301e+01_r8,      &
       5.2404e+02_r8, 2.3304e+03_r8, 1.0546e+04_r8, 2.5451e+03_r8,-3.2994e+04_r8,      &
      -1.8650e+03_r8, 4.9665e+03_r8, 1.8239e+03_r8, 3.6168e+01_r8, 6.6265e+02_r8,      &
       2.6543e+03_r8, 1.0920e+04_r8, 2.3793e+03_r8,-3.4430e+04_r8,-2.1422e+03_r8,      &
       5.3802e+03_r8, 1.8279e+03_r8, 4.0172e+01_r8, 8.1883e+02_r8, 2.9511e+03_r8,      &
       1.1315e+04_r8, 2.2543e+03_r8,-3.6037e+04_r8,-2.4455e+03_r8, 5.8820e+03_r8,      &
       1.8575e+03_r8, 4.4122e+01_r8, 9.9436e+02_r8, 3.1962e+03_r8, 1.1812e+04_r8,      &
       2.1571e+03_r8,-3.7838e+04_r8,-2.7947e+03_r8, 6.4488e+03_r8, 1.9130e+03_r8/

   data ((a360(ix,js),js=1,9),ix=16,20)/                                &
       4.8635e+01_r8, 1.1934e+03_r8, 3.3757e+03_r8, 1.2365e+04_r8, 2.0383e+03_r8,      &
      -3.9649e+04_r8,-3.1911e+03_r8, 7.0180e+03_r8, 1.9862e+03_r8, 5.4146e+01_r8,      &
       1.4193e+03_r8, 3.5325e+03_r8, 1.2990e+04_r8, 1.9202e+03_r8,-4.1577e+04_r8,      &
      -3.6210e+03_r8, 7.5891e+03_r8, 2.0823e+03_r8, 6.0556e+01_r8, 1.6705e+03_r8,      &
       3.6979e+03_r8, 1.3664e+04_r8, 1.8183e+03_r8,-4.3643e+04_r8,-4.0505e+03_r8,      &
       8.1520e+03_r8, 2.1988e+03_r8, 8.8450e+01_r8, 1.8763e+03_r8, 3.9281e+03_r8,      &
       1.4427e+04_r8, 1.7944e+03_r8,-4.6000e+04_r8,-4.4724e+03_r8, 8.7616e+03_r8,      &
       2.3413e+03_r8, 1.3066e+02_r8, 2.0515e+03_r8, 4.2101e+03_r8, 1.5254e+04_r8,      &
       1.8910e+03_r8,-4.8678e+04_r8,-4.8678e+03_r8, 9.4122e+03_r8, 2.5201e+03_r8/

   data ((a360(ix,js),js=1,9),ix=21,25)/                                &
       1.9627e+02_r8, 2.1888e+03_r8, 4.5373e+03_r8, 1.6134e+04_r8, 2.1484e+03_r8,      &
      -5.1718e+04_r8,-5.2303e+03_r8, 1.0080e+04_r8, 2.7544e+03_r8, 2.9253e+02_r8,      &
       2.2999e+03_r8, 4.9057e+03_r8, 1.7119e+04_r8, 2.5685e+03_r8,-5.5267e+04_r8,      &
      -5.5588e+03_r8, 1.0805e+04_r8, 3.0452e+03_r8, 4.2069e+02_r8, 2.4046e+03_r8,      &
       5.3221e+03_r8, 1.8287e+04_r8, 3.1724e+03_r8,-5.9480e+04_r8,-5.8650e+03_r8,      &
       1.1590e+04_r8, 3.4422e+03_r8, 5.6565e+02_r8, 2.5471e+03_r8, 5.7802e+03_r8,      &
       1.9750e+04_r8, 3.9225e+03_r8,-6.4671e+04_r8,-6.1721e+03_r8, 1.2437e+04_r8,      &
       4.0778e+03_r8, 6.8164e+02_r8, 2.7280e+03_r8, 6.2513e+03_r8, 2.1474e+04_r8,      &
       4.7161e+03_r8,-7.0534e+04_r8,-6.4551e+03_r8, 1.3250e+04_r8, 5.0017e+03_r8/

   data ((a360(ix,js),js=1,9),ix=26,31)/                                &
       7.5815e+02_r8, 2.9531e+03_r8, 6.7397e+03_r8, 2.3465e+04_r8, 5.5030e+03_r8,      &
      -7.7007e+04_r8,-6.7321e+03_r8, 1.4003e+04_r8, 6.2391e+03_r8, 8.1943e+02_r8,      &
       3.2232e+03_r8, 7.2806e+03_r8, 2.5738e+04_r8, 6.2592e+03_r8,-8.4190e+04_r8,      &
      -7.0625e+03_r8, 1.4724e+04_r8, 7.6781e+03_r8, 8.8343e+02_r8, 3.5281e+03_r8,      &
       7.9013e+03_r8, 2.8287e+04_r8, 6.9149e+03_r8,-9.2081e+04_r8,-7.4700e+03_r8,      &
       1.5391e+04_r8, 9.1934e+03_r8, 9.6259e+02_r8, 3.8750e+03_r8, 8.6310e+03_r8,      &
       3.0544e+04_r8, 7.4519e+03_r8,-9.9873e+04_r8,-7.8232e+03_r8, 1.5896e+04_r8,      &
       1.0647e+04_r8, 1.0657e+03_r8, 4.2669e+03_r8, 9.4391e+03_r8, 3.1950e+04_r8,      &
       7.2098e+03_r8,-1.0632e+05_r8,-7.8389e+03_r8, 1.6552e+04_r8, 1.1656e+04_r8,      &
       1.1918e+03_r8, 4.7413e+03_r8, 1.0192e+04_r8, 2.9994e+04_r8, 6.6001e+03_r8,      &
      -1.0795e+05_r8,-7.1220e+03_r8, 1.7021e+04_r8, 1.1715e+04_r8/

   data ((a360(ix,js),js=1,9),ix=32,37)/                                &
       1.3448e+03_r8, 5.3049e+03_r8, 1.0804e+04_r8, 2.8152e+04_r8, 3.0057e+03_r8,      &
      -1.0807e+05_r8,-5.8296e+03_r8, 1.9545e+04_r8, 1.0393e+04_r8, 1.5321e+03_r8,      &
       6.0282e+03_r8, 1.1182e+04_r8, 2.4855e+04_r8, 3.6060e+03_r8,-1.0871e+05_r8,      &
      -5.2856e+03_r8, 2.0777e+04_r8, 8.2516e+03_r8, 1.7158e+03_r8, 6.8949e+03_r8,      &
       1.0920e+04_r8, 2.3437e+04_r8, 1.3298e+03_r8,-1.1132e+05_r8,-4.5135e+03_r8,      &
       2.5203e+04_r8, 6.1970e+03_r8, 1.9016e+03_r8, 7.9667e+03_r8, 9.9618e+03_r8,      &
       2.0525e+04_r8, 4.2817e+03_r8,-1.1362e+05_r8,-5.0031e+03_r8, 2.7194e+04_r8,      &
       4.4984e+03_r8, 2.1244e+03_r8, 9.3828e+03_r8, 7.5590e+03_r8, 1.9678e+04_r8,      &
       7.0220e+03_r8,-1.2033e+05_r8,-4.9126e+03_r8, 3.1851e+04_r8, 3.8717e+03_r8,      &
       2.3664e+03_r8, 1.1265e+04_r8, 5.6693e+03_r8, 2.1491e+04_r8, 1.1265e+04_r8,      &
      -1.3697e+05_r8,-5.7812e+03_r8, 4.0486e+04_r8, 2.9647e+03_r8/

   data ((a360(ix,js),js=1,9),ix=38,43)/                                &
       2.7913e+03_r8, 1.4030e+04_r8, 4.8102e+03_r8, 2.6050e+04_r8, 1.3722e+04_r8,      &
      -1.7247e+05_r8,-5.1546e+03_r8, 5.9806e+04_r8, 4.0151e+03_r8, 3.2048e+03_r8,      &
       1.6855e+04_r8, 3.8414e+03_r8, 3.0267e+04_r8, 1.8679e+04_r8,-2.0901e+05_r8,      &
      -5.4977e+03_r8, 7.7718e+04_r8, 5.9062e+03_r8, 3.4514e+03_r8, 1.8244e+04_r8,      &
       1.8858e+03_r8, 3.4031e+04_r8, 2.2001e+04_r8,-2.3196e+05_r8,-6.1386e+03_r8,      &
       8.9356e+04_r8, 7.6331e+03_r8, 3.1735e+03_r8, 1.7858e+04_r8, 1.5239e+03_r8,      &
       3.9874e+04_r8, 3.0144e+04_r8,-2.3702e+05_r8,-1.2064e+04_r8, 8.4326e+04_r8,      &
       5.9467e+03_r8, 2.8899e+03_r8, 1.7226e+04_r8, 3.6062e+03_r8, 5.0991e+04_r8,      &
       3.8346e+04_r8,-2.5185e+05_r8,-1.9277e+04_r8, 8.1480e+04_r8, 4.3808e+03_r8,      &
       2.6308e+03_r8, 1.9510e+04_r8, 6.6321e+03_r8, 5.9251e+04_r8, 3.5217e+04_r8,      &
      -2.6413e+05_r8,-2.5896e+04_r8, 8.2704e+04_r8, 2.2345e+03_r8/

   data ((b360(ix,js),js=1,9),ix=1,5)/                                  &
       3.3220e+03_r8, 0.0000e+00_r8, 7.3567e+03_r8, 7.6803e+04_r8, 3.6454e+04_r8,      &
      -2.4080e+05_r8, 4.5492e+03_r8, 1.9812e+04_r8, 2.4305e+04_r8, 3.2356e+03_r8,      &
       0.0000e+00_r8, 8.9320e+03_r8, 7.6994e+04_r8, 3.5834e+04_r8,-2.4231e+05_r8,      &
       2.8598e+03_r8, 2.5713e+04_r8, 2.2445e+04_r8, 3.4175e+03_r8, 0.0000e+00_r8,      &
       1.2232e+04_r8, 7.7599e+04_r8, 2.9999e+04_r8,-2.4245e+05_r8,-8.9565e+02_r8,      &
       2.5185e+04_r8, 2.2348e+04_r8, 3.8227e+03_r8, 0.0000e+00_r8, 1.8057e+04_r8,      &
       7.7878e+04_r8, 2.1680e+04_r8,-2.4276e+05_r8,-4.7738e+03_r8, 2.1046e+04_r8,      &
       2.3454e+04_r8, 4.1639e+03_r8, 0.0000e+00_r8, 2.5061e+04_r8, 7.5590e+04_r8,      &
       1.6597e+04_r8,-2.4395e+05_r8,-7.3727e+03_r8, 1.8454e+04_r8, 2.4034e+04_r8/

   data ((b360(ix,js),js=1,9),ix=6,10)/                                 &
       2.0226e+03_r8, 3.6553e+03_r8, 2.7663e+04_r8, 7.4402e+04_r8, 1.3451e+04_r8,      &
      -2.4500e+05_r8,-9.0677e+03_r8, 1.7600e+04_r8, 2.4039e+04_r8, 2.3371e+03_r8,      &
       4.9460e+03_r8, 2.6742e+04_r8, 7.3298e+04_r8, 1.1572e+04_r8,-2.4282e+05_r8,      &
      -1.0360e+04_r8, 1.8247e+04_r8, 2.3128e+04_r8, 2.5821e+03_r8, 6.5452e+03_r8,      &
       2.3738e+04_r8, 7.1130e+04_r8, 9.2331e+03_r8,-2.3376e+05_r8,-1.1453e+04_r8,      &
       1.8905e+04_r8, 2.1477e+04_r8, 2.8407e+03_r8, 8.7467e+03_r8, 2.1169e+04_r8,      &
       6.9967e+04_r8, 6.3343e+03_r8,-2.2645e+05_r8,-1.2747e+04_r8, 1.9798e+04_r8,      &
       2.0046e+04_r8, 3.1777e+03_r8, 1.1663e+04_r8, 2.0303e+04_r8, 7.1744e+04_r8,      &
       3.4971e+03_r8,-2.2845e+05_r8,-1.4826e+04_r8, 2.1600e+04_r8, 1.9387e+04_r8/

   data ((b360(ix,js),js=1,9),ix=11,15)/                                &
       3.5824e+03_r8, 1.4759e+04_r8, 2.0870e+04_r8, 7.5297e+04_r8, 1.1064e+03_r8,      & 
      -2.3739e+05_r8,-1.7845e+04_r8, 2.4720e+04_r8, 1.8972e+04_r8, 4.0208e+03_r8,      &
       1.7477e+04_r8, 2.2103e+04_r8, 7.8013e+04_r8,-1.3374e+03_r8,-2.4578e+05_r8,      &
      -2.1347e+04_r8, 2.7998e+04_r8, 1.8416e+04_r8, 4.4998e+03_r8, 1.9891e+04_r8,      &
       2.3925e+04_r8, 8.0263e+04_r8,-4.0759e+03_r8,-2.5531e+05_r8,-2.4750e+04_r8,      &
       3.1873e+04_r8, 1.7813e+04_r8, 4.9537e+03_r8, 2.1855e+04_r8, 2.5057e+04_r8,      &
       8.0524e+04_r8,-6.6805e+03_r8,-2.5945e+05_r8,-2.7552e+04_r8, 3.5276e+04_r8,      &
       1.6945e+04_r8, 5.3966e+03_r8, 2.3717e+04_r8, 2.5315e+04_r8, 7.9668e+04_r8,      &
      -8.9218e+03_r8,-2.5948e+05_r8,-2.9806e+04_r8, 3.8224e+04_r8, 1.5980e+04_r8/

   data ((b360(ix,js),js=1,9),ix=16,20)/                                &
       5.8192e+03_r8, 2.5593e+04_r8, 2.5201e+04_r8, 8.0720e+04_r8,-9.4973e+03_r8,      &
      -2.6599e+05_r8,-3.1836e+04_r8, 4.2603e+04_r8, 1.5458e+04_r8, 6.2188e+03_r8,      &
       2.7630e+04_r8, 2.4953e+04_r8, 8.2660e+04_r8,-9.0242e+03_r8,-2.7596e+05_r8,      &
      -3.3635e+04_r8, 4.7478e+04_r8, 1.5519e+04_r8, 6.5500e+03_r8, 2.9757e+04_r8,      &
       2.5013e+04_r8, 8.6318e+04_r8,-7.6888e+03_r8,-2.9146e+05_r8,-3.5697e+04_r8,      &
       5.3080e+04_r8, 1.6177e+04_r8, 8.8577e+03_r8, 2.6042e+04_r8, 2.7715e+04_r8,      &
       9.1263e+04_r8,-5.6446e+03_r8,-3.1026e+05_r8,-3.8081e+04_r8, 5.8972e+04_r8,      &
       1.7414e+04_r8, 1.2119e+04_r8, 2.1474e+04_r8, 3.1224e+04_r8, 9.7933e+04_r8,      &
      -3.2719e+03_r8,-3.3360e+05_r8,-4.1185e+04_r8, 6.4976e+04_r8, 1.9369e+04_r8/

   data ((b360(ix,js),js=1,9),ix=21,25)/                                &
       1.6921e+04_r8, 1.6598e+04_r8, 3.5627e+04_r8, 1.0752e+05_r8,-5.1939e+02_r8,      &
      -3.6518e+05_r8,-4.5432e+04_r8, 7.1580e+04_r8, 2.2536e+04_r8, 2.3498e+04_r8,      &
       1.2022e+04_r8, 4.0689e+04_r8, 1.2032e+05_r8, 2.3753e+03_r8,-4.0492e+05_r8,      &
      -5.0818e+04_r8, 7.8740e+04_r8, 2.7036e+04_r8, 3.1012e+04_r8, 8.2263e+03_r8,      &
       4.6120e+04_r8, 1.3578e+05_r8, 5.0895e+03_r8,-4.4899e+05_r8,-5.7140e+04_r8,      &
       8.5635e+04_r8, 3.2388e+04_r8, 3.6715e+04_r8, 5.6657e+03_r8, 5.0954e+04_r8,      &
       1.5158e+05_r8, 6.0793e+03_r8,-4.8814e+05_r8,-6.2938e+04_r8, 8.8307e+04_r8,      &
       3.9590e+04_r8, 4.0345e+04_r8, 4.5077e+03_r8, 5.6889e+04_r8, 1.7253e+05_r8,      &
       4.4000e+03_r8,-5.3621e+05_r8,-7.0080e+04_r8, 8.9426e+04_r8, 5.0068e+04_r8/

   data ((b360(ix,js),js=1,9),ix=26,31)/                                &
       4.3147e+04_r8, 4.5945e+03_r8, 6.5637e+04_r8, 2.0208e+05_r8,-1.4353e+03_r8,      &
      -6.0197e+05_r8,-8.0122e+04_r8, 9.0303e+04_r8, 6.4638e+04_r8, 4.5803e+04_r8,      &
       6.1820e+03_r8, 7.8481e+04_r8, 2.4208e+05_r8,-1.2292e+04_r8,-6.9085e+05_r8,      &
      -9.4005e+04_r8, 9.2021e+04_r8, 8.2655e+04_r8, 4.8930e+04_r8, 9.7251e+03_r8,      &
       9.6788e+04_r8, 2.9438e+05_r8,-2.8855e+04_r8,-8.0949e+05_r8,-1.1269e+05_r8,      &
       9.5309e+04_r8, 1.0409e+05_r8, 5.2153e+04_r8, 1.5736e+04_r8, 1.2112e+05_r8,      &
       3.5047e+05_r8,-4.9630e+04_r8,-9.4867e+05_r8,-1.3405e+05_r8, 9.9646e+04_r8,      &
       1.2758e+05_r8, 5.5289e+04_r8, 2.5000e+04_r8, 1.5248e+05_r8, 4.0443e+05_r8,      &
      -7.8986e+04_r8,-1.1015e+06_r8,-1.5585e+05_r8, 1.0908e+05_r8, 1.4991e+05_r8,      &
       5.8421e+04_r8, 3.7738e+04_r8, 1.8933e+05_r8, 4.1350e+05_r8,-1.1087e+05_r8,      &
      -1.2104e+06_r8,-1.6912e+05_r8, 1.1441e+05_r8, 1.6453e+05_r8/

   data ((b360(ix,js),js=1,9),ix=32,37)/                                &
       6.2520e+04_r8, 5.6218e+04_r8, 2.3453e+05_r8, 4.2779e+05_r8,-1.7334e+05_r8,      &
      -1.3347e+06_r8,-1.8177e+05_r8, 1.4998e+05_r8, 1.6136e+05_r8, 6.7847e+04_r8,      &
       8.0624e+04_r8, 2.7957e+05_r8, 3.9240e+05_r8,-1.8747e+05_r8,-1.4391e+06_r8,      &
      -1.9884e+05_r8, 1.7445e+05_r8, 1.3890e+05_r8, 7.4315e+04_r8, 1.1285e+05_r8,      &   
       3.1347e+05_r8, 3.8243e+05_r8,-2.4294e+05_r8,-1.5620e+06_r8,-2.2047e+05_r8,      &
       2.3239e+05_r8, 1.1589e+05_r8, 8.1983e+04_r8, 1.5529e+05_r8, 3.2509e+05_r8,      &
       3.1830e+05_r8,-2.2382e+05_r8,-1.6485e+06_r8,-2.4782e+05_r8, 2.6696e+05_r8,      &
       9.0607e+04_r8, 9.1461e+04_r8, 2.1130e+05_r8, 2.8260e+05_r8, 2.8233e+05_r8,      &   
      -2.0994e+05_r8,-1.7623e+06_r8,-2.6761e+05_r8, 3.2415e+05_r8, 7.9589e+04_r8,      &
       1.0184e+05_r8, 2.8912e+05_r8, 2.4285e+05_r8, 2.8782e+05_r8,-1.9263e+05_r8,      &
      -2.0088e+06_r8,-3.0816e+05_r8, 4.3060e+05_r8, 6.2821e+04_r8/

   data ((b360(ix,js),js=1,9),ix=38,43)/                                &
       1.1806e+05_r8, 4.0223e+05_r8, 2.2517e+05_r8, 3.3973e+05_r8,-2.3078e+05_r8,      &
      -2.5372e+06_r8,-3.5716e+05_r8, 6.5288e+05_r8, 7.7642e+04_r8, 1.3586e+05_r8,      &
       5.2631e+05_r8, 1.8982e+05_r8, 3.9760e+05_r8,-2.1883e+05_r8,-3.1397e+06_r8,      &
      -4.0467e+05_r8, 8.8836e+05_r8, 1.0255e+05_r8, 1.5133e+05_r8, 6.1562e+05_r8,      &
       1.0783e+05_r8, 4.8241e+05_r8,-1.6362e+05_r8,-3.6844e+06_r8,-4.1372e+05_r8,      &
       1.1113e+06_r8, 1.2212e+05_r8, 1.5195e+05_r8, 6.5082e+05_r8, 7.2539e+04_r8,      &
       6.3688e+05_r8, 5.3880e+04_r8,-4.1600e+06_r8,-4.5347e+05_r8, 1.2093e+06_r8,      &
       9.2240e+04_r8, 1.4336e+05_r8, 6.1671e+05_r8, 1.2423e+05_r8, 8.1556e+05_r8,      &
       1.3919e+05_r8,-4.2989e+06_r8,-5.6767e+05_r8, 1.1435e+06_r8, 6.9843e+04_r8,      &
       1.3071e+05_r8, 6.5514e+05_r8, 1.9799e+05_r8, 8.8504e+05_r8, 9.0991e+04_r8,      &
      -4.2280e+06_r8,-6.1482e+05_r8, 1.0996e+06_r8, 3.0361e+04_r8/

   data ((a540(ix,js),js=1,9),ix=1,5)/                                  &
       1.4046e+01_r8, 0.0000e+00_r8, 8.4121e+01_r8, 2.7000e+03_r8, 4.9855e+03_r8,      &
      -1.4428e+04_r8, 1.1932e+03_r8, 3.1112e+03_r8, 5.2292e+02_r8, 1.8382e+01_r8,      &
       0.0000e+00_r8, 1.3808e+02_r8, 3.5508e+03_r8, 5.2498e+03_r8,-1.7036e+04_r8,      &
       1.1096e+03_r8, 3.6461e+03_r8, 7.1571e+02_r8, 2.2854e+01_r8, 0.0000e+00_r8,      &
       2.2670e+02_r8, 4.5327e+03_r8, 5.5670e+03_r8,-2.0072e+04_r8, 9.0698e+02_r8,      &
       4.0662e+03_r8, 9.6773e+02_r8, 2.6788e+01_r8, 0.0000e+00_r8, 3.6639e+02_r8,      &
       5.5619e+03_r8, 5.8004e+03_r8,-2.3126e+04_r8, 6.7679e+02_r8, 4.2755e+03_r8,      &
       1.2365e+03_r8, 3.1310e+01_r8, 0.0000e+00_r8, 5.8732e+02_r8, 6.6192e+03_r8,      &
       5.8620e+03_r8,-2.5990e+04_r8, 3.3715e+02_r8, 4.3206e+03_r8, 1.5084e+03_r8/

   data ((a540(ix,js),js=1,9),ix=6,10)/                                 &
       6.4043e+00_r8, 4.6835e+01_r8, 8.5364e+02_r8, 7.7821e+03_r8, 5.6476e+03_r8,      &
      -2.8674e+04_r8,-3.2399e+01_r8, 4.3542e+03_r8, 1.7554e+03_r8, 9.6264e+00_r8,      &
       7.8879e+01_r8, 1.1548e+03_r8, 8.8791e+03_r8, 5.2330e+03_r8,-3.1025e+04_r8,      &
      -3.5404e+02_r8, 4.4739e+03_r8, 1.9480e+03_r8, 1.3550e+01_r8, 1.3075e+02_r8,      &
       1.4391e+03_r8, 9.8097e+03_r8, 4.8032e+03_r8,-3.2995e+04_r8,-6.8661e+02_r8,      &
       4.6960e+03_r8, 2.0777e+03_r8, 1.7714e+01_r8, 2.0669e+02_r8, 1.6861e+03_r8,      &
       1.0641e+04_r8, 4.3793e+03_r8,-3.4751e+04_r8,-1.0292e+03_r8, 5.0045e+03_r8,      &
       2.1496e+03_r8, 2.1444e+01_r8, 3.0070e+02_r8, 1.9210e+03_r8, 1.1366e+04_r8,      &
       3.9995e+03_r8,-3.6335e+04_r8,-1.3571e+03_r8, 5.3448e+03_r8, 2.1701e+03_r8/

   data ((a540(ix,js),js=1,9),ix=11,15)/                                &
       2.4860e+01_r8, 4.0590e+02_r8, 2.1813e+03_r8, 1.1988e+04_r8, 3.7302e+03_r8,      &
      -3.7857e+04_r8,-1.6700e+03_r8, 5.7183e+03_r8, 2.1533e+03_r8, 2.8308e+01_r8,      &
       5.2176e+02_r8, 2.5305e+03_r8, 1.2550e+04_r8, 3.5359e+03_r8,-3.9488e+04_r8,      &
      -2.0723e+03_r8, 6.1238e+03_r8, 2.1332e+03_r8, 3.2231e+01_r8, 6.6569e+02_r8,      &
       2.9583e+03_r8, 1.3080e+04_r8, 3.2613e+03_r8,-4.1264e+04_r8,-2.4406e+03_r8,      &
       6.6123e+03_r8, 2.1193e+03_r8, 3.6139e+01_r8, 8.3147e+02_r8, 3.3614e+03_r8,      &
       1.3578e+04_r8, 3.0508e+03_r8,-4.3167e+04_r8,-2.8085e+03_r8, 7.2030e+03_r8,      &
       2.1301e+03_r8, 4.0256e+01_r8, 1.0297e+03_r8, 3.7022e+03_r8, 1.4171e+04_r8,      &
       2.9128e+03_r8,-4.5336e+04_r8,-3.2104e+03_r8, 7.8946e+03_r8, 2.1776e+03_r8/

   data ((a540(ix,js),js=1,9),ix=16,20)/                                &
       4.4889e+01_r8, 1.2599e+03_r8, 3.9615e+03_r8, 1.4830e+04_r8, 2.8030e+03_r8,      &
      -4.7625e+04_r8,-3.6641e+03_r8, 8.6225e+03_r8, 2.2558e+03_r8, 5.0015e+01_r8,      &
       1.5224e+03_r8, 4.1757e+03_r8, 1.5557e+04_r8, 2.7204e+03_r8,-5.0044e+04_r8,      &
      -4.1690e+03_r8, 9.3569e+03_r8, 2.3618e+03_r8, 5.6369e+01_r8, 1.8214e+03_r8,      &
       4.3850e+03_r8, 1.6354e+04_r8, 2.6560e+03_r8,-5.2643e+04_r8,-4.6917e+03_r8,      &
       1.0073e+04_r8, 2.4979e+03_r8, 8.3103e+01_r8, 2.0790e+03_r8, 4.6427e+03_r8,      &
       1.7254e+04_r8, 2.6678e+03_r8,-5.5501e+04_r8,-5.2101e+03_r8, 1.0809e+04_r8,      &
       2.6617e+03_r8, 1.2393e+02_r8, 2.3045e+03_r8, 4.9515e+03_r8, 1.8262e+04_r8,      &
       2.8040e+03_r8,-5.8738e+04_r8,-5.7019e+03_r8, 1.1579e+04_r8, 2.8680e+03_r8/

   data ((a540(ix,js),js=1,9),ix=21,25)/                                &
       1.8908e+02_r8, 2.4870e+03_r8, 5.3135e+03_r8, 1.9357e+04_r8, 3.1110e+03_r8,      &
      -6.2419e+04_r8,-6.1562e+03_r8, 1.2371e+04_r8, 3.1342e+03_r8, 2.8397e+02_r8,      &
       2.6277e+03_r8, 5.7304e+03_r8, 2.0564e+04_r8, 3.6012e+03_r8,-6.6647e+04_r8,      &
      -6.5563e+03_r8, 1.3212e+04_r8, 3.4531e+03_r8, 4.1265e+02_r8, 2.7604e+03_r8,      &   
       6.2057e+03_r8, 2.1959e+04_r8, 4.3164e+03_r8,-7.1623e+04_r8,-6.9163e+03_r8,      &
       1.4129e+04_r8, 3.8711e+03_r8, 5.5908e+02_r8, 2.9304e+03_r8, 6.7324e+03_r8,      &
       2.3661e+04_r8, 5.2276e+03_r8,-7.7697e+04_r8,-7.2578e+03_r8, 1.5148e+04_r8,      &
       4.5144e+03_r8, 6.7735e+02_r8, 3.1384e+03_r8, 7.2790e+03_r8, 2.5647e+04_r8,      &
       6.2256e+03_r8,-8.4578e+04_r8,-7.5573e+03_r8, 1.6206e+04_r8, 5.4499e+03_r8/

   data ((a540(ix,js),js=1,9),ix=26,31)/                                &
       7.6095e+02_r8, 3.3936e+03_r8, 7.8572e+03_r8, 2.7951e+04_r8, 7.2729e+03_r8,      &
      -9.2312e+04_r8,-7.8460e+03_r8, 1.7305e+04_r8, 6.7468e+03_r8, 8.2152e+02_r8,      &
       3.6897e+03_r8, 8.4910e+03_r8, 3.0565e+04_r8, 8.3712e+03_r8,-1.0096e+05_r8,      &
      -8.1788e+03_r8, 1.8426e+04_r8, 8.3335e+03_r8, 8.9734e+02_r8, 4.0371e+03_r8,      &
       9.2179e+03_r8, 3.3541e+04_r8, 9.4243e+03_r8,-1.1068e+05_r8,-8.6295e+03_r8,      &
       1.9539e+04_r8, 1.0194e+04_r8, 9.8654e+02_r8, 4.4310e+03_r8, 1.0051e+04_r8,      &
       3.6202e+04_r8, 1.0341e+04_r8,-1.2033e+05_r8,-9.0465e+03_r8, 2.0323e+04_r8,      &
       1.2170e+04_r8, 1.1010e+03_r8, 4.8898e+03_r8, 1.0979e+04_r8, 3.8051e+04_r8,      &
       1.0183e+04_r8,-1.2865e+05_r8,-9.1689e+03_r8, 2.1141e+04_r8, 1.3739e+04_r8,      &
       1.2523e+03_r8, 5.4507e+03_r8, 1.1845e+04_r8, 3.6056e+04_r8, 9.2641e+03_r8,      &
      -1.3103e+05_r8,-8.5702e+03_r8, 2.1465e+04_r8, 1.4032e+04_r8/

   data ((a540(ix,js),js=1,9),ix=32,37)/                                &
       1.4290e+03_r8, 6.1493e+03_r8, 1.2630e+04_r8, 3.4660e+04_r8, 4.6113e+03_r8,      &
      -1.3253e+05_r8,-7.6923e+03_r8, 2.5385e+04_r8, 1.1814e+04_r8, 1.6393e+03_r8,      &
       6.9603e+03_r8, 1.3045e+04_r8, 3.1292e+04_r8, 5.0250e+03_r8,-1.3241e+05_r8,      &
      -8.0218e+03_r8, 2.7248e+04_r8, 8.1799e+03_r8, 1.8422e+03_r8, 7.9125e+03_r8,      &
       1.2850e+04_r8, 2.9919e+04_r8, 1.6539e+03_r8,-1.3374e+05_r8,-7.5523e+03_r8,      &
       3.2067e+04_r8, 5.4650e+03_r8, 2.0556e+03_r8, 9.1198e+03_r8, 1.1952e+04_r8,      &
       2.6229e+04_r8, 5.4064e+03_r8,-1.3513e+05_r8,-8.5919e+03_r8, 3.3805e+04_r8,      &
       3.3909e+03_r8, 2.3195e+03_r8, 1.0765e+04_r8, 9.4232e+03_r8, 2.4680e+04_r8,      &
       7.8282e+03_r8,-1.4260e+05_r8,-8.1716e+03_r8, 3.9314e+04_r8, 3.2746e+03_r8,      &
       2.6741e+03_r8, 1.3056e+04_r8, 7.3683e+03_r8, 2.5901e+04_r8, 1.1853e+04_r8,      &
      -1.6172e+05_r8,-8.6679e+03_r8, 4.9358e+04_r8, 3.1283e+03_r8/

   data ((a540(ix,js),js=1,9),ix=38,43)/                                &
       3.1274e+03_r8, 1.6190e+04_r8, 6.6171e+03_r8, 2.9852e+04_r8, 1.3996e+04_r8,      &
      -2.0033e+05_r8,-7.1656e+03_r8, 7.0586e+04_r8, 5.0072e+03_r8, 3.5342e+03_r8,      &
       1.9357e+04_r8, 5.6559e+03_r8, 3.3611e+04_r8, 1.7828e+04_r8,-2.3702e+05_r8,      &
      -6.2284e+03_r8, 8.8163e+04_r8, 7.3170e+03_r8, 3.7834e+03_r8, 2.0965e+04_r8,      &
       3.2588e+03_r8, 3.7010e+04_r8, 2.1189e+04_r8,-2.5903e+05_r8,-5.8621e+03_r8,      &
       9.8554e+04_r8, 8.6223e+03_r8, 3.5046e+03_r8, 2.0797e+04_r8, 1.9896e+03_r8,      &
       4.3329e+04_r8, 3.0815e+04_r8,-2.6553e+05_r8,-1.1030e+04_r8, 9.3433e+04_r8,      &
       6.0978e+03_r8, 3.3389e+03_r8, 2.0607e+04_r8, 2.5567e+03_r8, 5.6928e+04_r8,      &
       4.3500e+04_r8,-2.8735e+05_r8,-1.7834e+04_r8, 8.9701e+04_r8, 3.7020e+03_r8,      &
       3.2675e+03_r8, 2.3717e+04_r8, 3.5488e+03_r8, 7.1174e+04_r8, 4.9418e+04_r8,      &
      -3.1866e+05_r8,-2.5099e+04_r8, 9.4302e+04_r8, 1.5043e+03_r8/

   data ((b540(ix,js),js=1,9),ix=1,5)/                                  &
       2.9709e+03_r8, 0.0000e+00_r8, 7.2631e+03_r8, 8.2400e+04_r8, 4.6777e+04_r8,      &
      -2.6896e+05_r8, 8.8736e+03_r8, 2.2636e+04_r8, 2.6978e+04_r8, 3.1238e+03_r8,      &
       0.0000e+00_r8, 9.6559e+03_r8, 8.7263e+04_r8, 4.5581e+04_r8,-2.8113e+05_r8,      &
       6.5338e+03_r8, 3.0570e+04_r8, 2.5612e+04_r8, 3.4594e+03_r8, 0.0000e+00_r8,      &
       1.3781e+04_r8, 9.1771e+04_r8, 4.2054e+04_r8,-2.9581e+05_r8, 2.4720e+03_r8,      &
       3.2871e+04_r8, 2.6461e+04_r8, 3.9369e+03_r8, 0.0000e+00_r8, 2.0602e+04_r8,      &
       9.1745e+04_r8, 3.0580e+04_r8,-2.9239e+05_r8,-2.5695e+03_r8, 2.6492e+04_r8,      &
       2.8006e+04_r8, 4.3495e+03_r8, 0.0000e+00_r8, 2.8814e+04_r8, 8.7356e+04_r8,      &
       2.1482e+04_r8,-2.8467e+05_r8,-6.3525e+03_r8, 2.1028e+04_r8, 2.8623e+04_r8/

   data ((b540(ix,js),js=1,9),ix=6,10)/                                 &
       2.0005e+03_r8, 4.0531e+03_r8, 3.1958e+04_r8, 8.4074e+04_r8, 1.6201e+04_r8,      &
      -2.7909e+05_r8,-8.6225e+03_r8, 1.8720e+04_r8, 2.8409e+04_r8, 2.3629e+03_r8,      &
       5.5993e+03_r8, 3.1028e+04_r8, 8.2417e+04_r8, 1.3922e+04_r8,-2.7599e+05_r8,      &
      -1.0282e+04_r8, 1.9284e+04_r8, 2.7425e+04_r8, 2.6727e+03_r8, 7.5414e+03_r8,      &
       2.7733e+04_r8, 8.1705e+04_r8, 1.1962e+04_r8,-2.7148e+05_r8,-1.1901e+04_r8,      &
       2.0897e+04_r8, 2.5922e+04_r8, 2.9924e+03_r8, 1.0195e+04_r8, 2.4967e+04_r8,      &
       8.2254e+04_r8, 9.0794e+03_r8,-2.6878e+05_r8,-1.3758e+04_r8, 2.2862e+04_r8,      &
       2.4609e+04_r8, 3.3628e+03_r8, 1.3635e+04_r8, 2.4057e+04_r8, 8.5280e+04_r8,      &
       5.7370e+03_r8,-2.7339e+05_r8,-1.6337e+04_r8, 2.5177e+04_r8, 2.4060e+04_r8/

   data ((b540(ix,js),js=1,9),ix=11,15)/                                &
       3.7910e+03_r8, 1.7353e+04_r8, 2.4740e+04_r8, 8.9874e+04_r8, 2.8386e+03_r8,      &
      -2.8470e+05_r8,-2.0176e+04_r8, 2.9246e+04_r8, 2.3552e+04_r8, 4.2537e+03_r8,      &
       2.0528e+04_r8, 2.6524e+04_r8, 9.4192e+04_r8,-1.9538e+02_r8,-2.9713e+05_r8,      &
      -2.4840e+04_r8, 3.3537e+04_r8, 2.3038e+04_r8, 4.7553e+03_r8, 2.3449e+04_r8,      &
       2.8849e+04_r8, 9.7935e+04_r8,-3.3681e+03_r8,-3.1112e+05_r8,-2.9373e+04_r8,      &
       3.8841e+04_r8, 2.2236e+04_r8, 5.2435e+03_r8, 2.5812e+04_r8, 3.0777e+04_r8,      &
       1.0004e+05_r8,-6.4184e+03_r8,-3.2110e+05_r8,-3.3387e+04_r8, 4.3876e+04_r8,      &
       2.1211e+04_r8, 5.7016e+03_r8, 2.8053e+04_r8, 3.1613e+04_r8, 1.0025e+05_r8,      &
      -9.2545e+03_r8,-3.2451e+05_r8,-3.6675e+04_r8, 4.8014e+04_r8, 1.9899e+04_r8/

   data ((b540(ix,js),js=1,9),ix=16,20)/                                &
       6.1442e+03_r8, 3.0360e+04_r8, 3.1866e+04_r8, 1.0151e+05_r8,-1.0762e+04_r8,      &
      -3.3154e+05_r8,-3.9524e+04_r8, 5.3002e+04_r8, 1.8817e+04_r8, 6.5885e+03_r8,      &
       3.3054e+04_r8, 3.1916e+04_r8, 1.0396e+05_r8,-1.0604e+04_r8,-3.4429e+05_r8,      &
      -4.1917e+04_r8, 5.8965e+04_r8, 1.8358e+04_r8, 6.9707e+03_r8, 3.5930e+04_r8,      &
       3.2193e+04_r8, 1.0817e+05_r8,-8.9122e+03_r8,-3.6374e+05_r8,-4.4279e+04_r8,      &
       6.5938e+04_r8, 1.8630e+04_r8, 9.5116e+03_r8, 3.2665e+04_r8, 3.5521e+04_r8,      &
       1.1452e+05_r8,-5.9078e+03_r8,-3.9000e+05_r8,-4.7090e+04_r8, 7.3913e+04_r8,      &
       1.9693e+04_r8, 1.3035e+04_r8, 2.8272e+04_r8, 3.9666e+04_r8, 1.2272e+05_r8,      &
      -1.9538e+03_r8,-4.2149e+05_r8,-5.0485e+04_r8, 8.2204e+04_r8, 2.1617e+04_r8/

   data ((b540(ix,js),js=1,9),ix=21,25)/                                &
       1.8173e+04_r8, 2.3394e+04_r8, 4.4584e+04_r8, 1.3380e+05_r8, 2.9300e+03_r8,      &
      -4.6142e+05_r8,-5.4863e+04_r8, 9.1181e+04_r8, 2.4808e+04_r8, 2.5309e+04_r8,      &
       1.8719e+04_r8, 5.0230e+04_r8, 1.4885e+05_r8, 8.6209e+03_r8,-5.1236e+05_r8,      &
      -6.0511e+04_r8, 1.0139e+05_r8, 2.9436e+04_r8, 3.3629e+04_r8, 1.4750e+04_r8,      &
       5.6231e+04_r8, 1.6736e+05_r8, 1.4927e+04_r8,-5.7045e+05_r8,-6.7277e+04_r8,      &
       1.1217e+05_r8, 3.5060e+04_r8, 3.9981e+04_r8, 1.1960e+04_r8, 6.1755e+04_r8,      &
       1.8744e+05_r8, 2.0137e+04_r8,-6.2674e+05_r8,-7.3834e+04_r8, 1.1934e+05_r8,      &
       4.3078e+04_r8, 4.4113e+04_r8, 1.0697e+04_r8, 6.8456e+04_r8, 2.1400e+05_r8,      &
       2.3263e+04_r8,-6.9595e+05_r8,-8.2136e+04_r8, 1.2566e+05_r8, 5.5471e+04_r8/

   data ((b540(ix,js),js=1,9),ix=26,31)/                                &
       4.7268e+04_r8, 1.0742e+04_r8, 7.7893e+04_r8, 2.5024e+05_r8, 2.2431e+04_r8,      &
      -7.8521e+05_r8,-9.3714e+04_r8, 1.3164e+05_r8, 7.3653e+04_r8, 5.0261e+04_r8,      &
       1.2320e+04_r8, 9.1514e+04_r8, 2.9871e+05_r8, 1.5721e+04_r8,-9.0003e+05_r8,      &
      -1.0984e+05_r8, 1.3732e+05_r8, 9.7259e+04_r8, 5.3382e+04_r8, 1.5997e+04_r8,      &
       1.1025e+05_r8, 3.6030e+05_r8, 1.6491e+03_r8,-1.0428e+06_r8,-1.3104e+05_r8,      &
       1.4263e+05_r8, 1.2572e+05_r8, 5.6536e+04_r8, 2.2075e+04_r8, 1.3502e+05_r8,      &
       4.2540e+05_r8,-1.9356e+04_r8,-1.2017e+06_r8,-1.5497e+05_r8, 1.4620e+05_r8,      &
       1.5667e+05_r8, 5.9668e+04_r8, 3.1332e+04_r8, 1.6712e+05_r8, 4.8760e+05_r8,      &
      -5.2497e+04_r8,-1.3721e+06_r8,-1.7921e+05_r8, 1.5491e+05_r8, 1.8500e+05_r8,      &
       6.3128e+04_r8, 4.4306e+04_r8, 2.0538e+05_r8, 5.0036e+05_r8,-8.8315e+04_r8,      &
      -1.5009e+06_r8,-1.9544e+05_r8, 1.6286e+05_r8, 2.0060e+05_r8/

   data ((b540(ix,js),js=1,9),ix=32,37)/                                &
       6.7252e+04_r8, 6.1350e+04_r8, 2.4981e+05_r8, 5.1866e+05_r8,-1.7393e+05_r8,      &
      -1.6079e+06_r8,-2.1616e+05_r8, 2.0514e+05_r8, 1.8421e+05_r8, 7.2523e+04_r8,      &
       8.7077e+04_r8, 3.0032e+05_r8, 4.9795e+05_r8,-1.9374e+05_r8,-1.7488e+06_r8,      &
      -2.5154e+05_r8, 2.5270e+05_r8, 1.4108e+05_r8, 7.8557e+04_r8, 1.2085e+05_r8,      &
       3.4369e+05_r8, 5.0559e+05_r8,-2.7072e+05_r8,-1.8995e+06_r8,-2.8633e+05_r8,      &
       3.2697e+05_r8, 1.0954e+05_r8, 8.6092e+04_r8, 1.6790e+05_r8, 3.6976e+05_r8,      &
       4.3787e+05_r8,-2.5173e+05_r8,-2.0235e+06_r8,-3.3252e+05_r8, 3.7066e+05_r8,      &
       7.9724e+04_r8, 9.6434e+04_r8, 2.3315e+05_r8, 3.3902e+05_r8, 3.9461e+05_r8,      &
      -2.5455e+05_r8,-2.2059e+06_r8,-3.6308e+05_r8, 4.4700e+05_r8, 8.0952e+04_r8,      &
       1.0946e+05_r8, 3.2289e+05_r8, 3.0583e+05_r8, 3.8625e+05_r8,-2.5579e+05_r8,      &
      -2.5156e+06_r8,-4.1451e+05_r8, 5.7253e+05_r8, 7.8418e+04_r8/

   data ((b540(ix,js),js=1,9),ix=38,43)/                                &
       1.2692e+05_r8, 4.4886e+05_r8, 3.0244e+05_r8, 4.2615e+05_r8,-3.1844e+05_r8,      &
      -3.1255e+06_r8,-4.6924e+05_r8, 8.2340e+05_r8, 1.1046e+05_r8, 1.4441e+05_r8,      &
       5.8243e+05_r8, 2.7541e+05_r8, 4.7182e+05_r8,-3.3649e+05_r8,-3.7185e+06_r8,      &
      -5.1067e+05_r8, 1.0440e+06_r8, 1.4277e+05_r8, 1.5916e+05_r8, 6.7402e+05_r8,      &
       1.8070e+05_r8, 5.5047e+05_r8,-2.9324e+05_r8,-4.2164e+06_r8,-5.1434e+05_r8,      &
       1.2427e+06_r8, 1.5864e+05_r8, 1.6012e+05_r8, 7.2380e+05_r8, 1.1475e+05_r8,      &
       7.1511e+05_r8,-6.9798e+04_r8,-4.6914e+06_r8,-5.4151e+05_r8, 1.3302e+06_r8,      &
       1.1012e+05_r8, 1.5471e+05_r8, 7.0419e+05_r8, 1.1310e+05_r8, 8.8568e+05_r8,      &
       2.4896e+04_r8,-4.6768e+06_r8,-6.3312e+05_r8, 1.1713e+06_r8, 6.4977e+04_r8,      &
       1.5021e+05_r8, 7.9025e+05_r8, 1.1779e+05_r8, 1.0507e+06_r8, 1.1771e+05_r8,      &
      -4.9507e+06_r8,-6.8570e+05_r8, 1.1955e+06_r8, 2.3632e+04_r8/

   data ((a720(ix,js),js=1,9),ix=1,5)/                                  &
       1.1000e+01_r8, 0.0000e+00_r8, 8.1568e+01_r8, 2.8044e+03_r8, 5.3898e+03_r8,      &
      -1.5286e+04_r8, 1.3891e+03_r8, 3.3052e+03_r8, 5.0976e+02_r8, 1.5350e+01_r8,      &
       0.0000e+00_r8, 1.3680e+02_r8, 3.6955e+03_r8, 5.7560e+03_r8,-1.8147e+04_r8,      &
       1.3596e+03_r8, 3.9557e+03_r8, 6.9200e+02_r8, 1.9863e+01_r8, 0.0000e+00_r8,      &
       2.2561e+02_r8, 4.7049e+03_r8, 6.2084e+03_r8,-2.1486e+04_r8, 1.1961e+03_r8,      &
       4.5137e+03_r8, 9.4594e+02_r8, 2.4344e+01_r8, 0.0000e+00_r8, 3.6564e+02_r8,      &
       5.7723e+03_r8, 6.6036e+03_r8,-2.4949e+04_r8, 9.9205e+02_r8, 4.8226e+03_r8,      &
       1.2417e+03_r8, 2.8832e+01_r8, 0.0000e+00_r8, 5.8203e+02_r8, 6.9163e+03_r8,      &
       6.8325e+03_r8,-2.8322e+04_r8, 6.4255e+02_r8, 4.9144e+03_r8, 1.5622e+03_r8/

   data ((a720(ix,js),js=1,9),ix=6,10)/                                 &
       4.8883e+00_r8, 4.4839e+01_r8, 8.4544e+02_r8, 8.2527e+03_r8, 6.7546e+03_r8,      &
      -3.1641e+04_r8, 2.3650e+02_r8, 4.9722e+03_r8, 1.8764e+03_r8, 7.5913e+00_r8,      &
       7.6632e+01_r8, 1.1494e+03_r8, 9.5596e+03_r8, 6.3700e+03_r8,-3.4538e+04_r8,      &
      -1.2944e+02_r8, 5.1143e+03_r8, 2.1263e+03_r8, 1.1092e+01_r8, 1.2916e+02_r8,      &
       1.4465e+03_r8, 1.0700e+04_r8, 5.9160e+03_r8,-3.6968e+04_r8,-5.1496e+02_r8,      &
       5.3844e+03_r8, 2.2995e+03_r8, 1.5051e+01_r8, 2.0677e+02_r8, 1.7088e+03_r8,      &
       1.1710e+04_r8, 5.4364e+03_r8,-3.9056e+04_r8,-9.2065e+02_r8, 5.7632e+03_r8,      &
       2.3914e+03_r8, 1.8844e+01_r8, 3.0261e+02_r8, 1.9641e+03_r8, 1.2606e+04_r8,      &
       4.9993e+03_r8,-4.0960e+04_r8,-1.3152e+03_r8, 6.1798e+03_r8, 2.4231e+03_r8/

   data ((a720(ix,js),js=1,9),ix=11,15)/                                &
       2.2391e+01_r8, 4.0801e+02_r8, 2.2543e+03_r8, 1.3383e+04_r8, 4.6738e+03_r8,      &
      -4.2763e+04_r8,-1.6992e+03_r8, 6.6211e+03_r8, 2.4071e+03_r8, 2.5971e+01_r8,      &
       5.2138e+02_r8, 2.6560e+03_r8, 1.4114e+04_r8, 4.4287e+03_r8,-4.4726e+04_r8,      &
      -2.1995e+03_r8, 7.0889e+03_r8, 2.3855e+03_r8, 2.9763e+01_r8, 6.6460e+02_r8,      &
       3.1664e+03_r8, 1.4809e+04_r8, 4.0645e+03_r8,-4.6843e+04_r8,-2.6537e+03_r8,      &
       7.6437e+03_r8, 2.3613e+03_r8, 3.3666e+01_r8, 8.3830e+02_r8, 3.6616e+03_r8,      &
       1.5430e+04_r8, 3.7727e+03_r8,-4.9054e+04_r8,-3.0909e+03_r8, 8.3082e+03_r8,      &
       2.3574e+03_r8, 3.7748e+01_r8, 1.0512e+03_r8, 4.0875e+03_r8, 1.6120e+04_r8,      &
       3.5818e+03_r8,-5.1529e+04_r8,-3.5443e+03_r8, 9.0933e+03_r8, 2.3920e+03_r8/

   data ((a720(ix,js),js=1,9),ix=16,20)/                                &
       4.2296e+01_r8, 1.3036e+03_r8, 4.4196e+03_r8, 1.6881e+04_r8, 3.4671e+03_r8,      &
      -5.4229e+04_r8,-4.0449e+03_r8, 9.9546e+03_r8, 2.4679e+03_r8, 4.7480e+01_r8,      &
       1.5965e+03_r8, 4.6912e+03_r8, 1.7692e+04_r8, 3.4059e+03_r8,-5.7053e+04_r8,      &
      -4.6035e+03_r8, 1.0823e+04_r8, 2.5794e+03_r8, 5.3627e+01_r8, 1.9302e+03_r8,      &
       4.9402e+03_r8, 1.8564e+04_r8, 3.3780e+03_r8,-6.0042e+04_r8,-5.1911e+03_r8,      &
       1.1665e+04_r8, 2.7244e+03_r8, 7.9167e+01_r8, 2.2283e+03_r8, 5.2269e+03_r8,      &
       1.9564e+04_r8, 3.4343e+03_r8,-6.3324e+04_r8,-5.7858e+03_r8, 1.2522e+04_r8,      &
       2.9022e+03_r8, 1.1948e+02_r8, 2.4959e+03_r8, 5.5584e+03_r8, 2.0703e+04_r8,      &
       3.6185e+03_r8,-6.7027e+04_r8,-6.3585e+03_r8, 1.3402e+04_r8, 3.1300e+03_r8/

   data ((a720(ix,js),js=1,9),ix=21,25)/                                &
       1.8391e+02_r8, 2.7146e+03_r8, 5.9441e+03_r8, 2.1967e+04_r8, 3.9800e+03_r8,      &
      -7.1233e+04_r8,-6.8930e+03_r8, 1.4300e+04_r8, 3.4242e+03_r8, 2.7904e+02_r8,      &
       2.8863e+03_r8, 6.3937e+03_r8, 2.3372e+04_r8, 4.5327e+03_r8,-7.6066e+04_r8,      &
      -7.3642e+03_r8, 1.5250e+04_r8, 3.7769e+03_r8, 4.0911e+02_r8, 3.0443e+03_r8,      &
       6.9146e+03_r8, 2.4981e+04_r8, 5.3338e+03_r8,-8.1711e+04_r8,-7.7825e+03_r8,      &
       1.6272e+04_r8, 4.2287e+03_r8, 5.5647e+02_r8, 3.2369e+03_r8, 7.4969e+03_r8,      &
       2.6902e+04_r8, 6.3584e+03_r8,-8.8522e+04_r8,-8.1632e+03_r8, 1.7400e+04_r8,      &
       4.8989e+03_r8, 6.7746e+02_r8, 3.4639e+03_r8, 8.1044e+03_r8, 2.9116e+04_r8,      &
       7.4937e+03_r8,-9.6201e+04_r8,-8.4797e+03_r8, 1.8601e+04_r8, 5.8474e+03_r8/

   data ((a720(ix,js),js=1,9),ix=26,31)/                                &
       7.6419e+02_r8, 3.7351e+03_r8, 8.7528e+03_r8, 3.1669e+04_r8, 8.7109e+03_r8,      &
      -1.0484e+05_r8,-8.7688e+03_r8, 1.9909e+04_r8, 7.1538e+03_r8, 8.3217e+02_r8,      &
       4.0555e+03_r8, 9.4681e+03_r8, 3.4560e+04_r8, 1.0047e+04_r8,-1.1462e+05_r8,      &
      -9.0931e+03_r8, 2.1351e+04_r8, 8.7935e+03_r8, 9.1301e+02_r8, 4.4297e+03_r8,      &
       1.0290e+04_r8, 3.7864e+04_r8, 1.1445e+04_r8,-1.2584e+05_r8,-9.5391e+03_r8,      &
       2.2897e+04_r8, 1.0786e+04_r8, 1.0164e+03_r8, 4.8721e+03_r8, 1.1249e+04_r8,      &
       4.0878e+04_r8, 1.2767e+04_r8,-1.3731e+05_r8,-9.9925e+03_r8, 2.4149e+04_r8,      &
       1.3084e+04_r8, 1.1415e+03_r8, 5.3815e+03_r8, 1.2293e+04_r8, 4.2941e+04_r8,      &
       1.2877e+04_r8,-1.4705e+05_r8,-1.0110e+04_r8, 2.5288e+04_r8, 1.5084e+04_r8,      &
       1.2975e+03_r8, 5.9819e+03_r8, 1.3226e+04_r8, 4.0585e+04_r8, 1.1739e+04_r8,      &
      -1.4918e+05_r8,-9.3569e+03_r8, 2.5526e+04_r8, 1.5795e+04_r8/

   data ((a720(ix,js),js=1,9),ix=32,37)/                                &
       1.4828e+03_r8, 6.7270e+03_r8, 1.4083e+04_r8, 3.9252e+04_r8, 6.2124e+03_r8,      &
      -1.5042e+05_r8,-8.4580e+03_r8, 2.9972e+04_r8, 1.3334e+04_r8, 1.7053e+03_r8,      &
       7.6200e+03_r8, 1.4682e+04_r8, 3.6183e+04_r8, 5.9149e+03_r8,-1.5116e+05_r8,      &
      -9.3969e+03_r8, 3.2817e+04_r8, 9.0065e+03_r8, 1.9276e+03_r8, 8.6904e+03_r8,      &
       1.4523e+04_r8, 3.5751e+04_r8, 2.4091e+03_r8,-1.5492e+05_r8,-1.0156e+04_r8,      &
       3.9419e+04_r8, 4.9338e+03_r8, 2.1451e+03_r8, 1.0040e+04_r8, 1.3698e+04_r8,      &
       3.1820e+04_r8, 7.1113e+03_r8,-1.5760e+05_r8,-1.1397e+04_r8, 3.9737e+04_r8,      &
       3.4762e+03_r8, 2.4469e+03_r8, 1.1859e+04_r8, 1.0944e+04_r8, 2.9494e+04_r8,      &
       8.8440e+03_r8,-1.6440e+05_r8,-9.9882e+03_r8, 4.4020e+04_r8, 4.6728e+03_r8,      &
       2.8411e+03_r8, 1.4376e+04_r8, 8.8420e+03_r8, 3.0261e+04_r8, 1.2589e+04_r8,      &
      -1.8408e+05_r8,-9.3312e+03_r8, 5.2682e+04_r8, 6.4506e+03_r8/

   data ((a720(ix,js),js=1,9),ix=38,43)/                                &
       3.3274e+03_r8, 1.7805e+04_r8, 8.3315e+03_r8, 3.4132e+04_r8, 1.3953e+04_r8,      &
      -2.2590e+05_r8,-7.2006e+03_r8, 7.6708e+04_r8, 8.4541e+03_r8, 3.7593e+03_r8,      &
       2.1296e+04_r8, 7.4382e+03_r8, 3.7788e+04_r8, 1.7525e+04_r8,-2.6251e+05_r8,      &
      -6.7751e+03_r8, 9.5562e+04_r8, 9.7409e+03_r8, 4.0215e+03_r8, 2.3146e+04_r8,      &
       4.4673e+03_r8, 3.9980e+04_r8, 1.9665e+04_r8,-2.8011e+05_r8,-6.2491e+03_r8,      &
       1.0617e+05_r8, 9.7124e+03_r8, 3.7615e+03_r8, 2.3278e+04_r8, 2.4201e+03_r8,      &
       4.5810e+04_r8, 2.9998e+04_r8,-2.8413e+05_r8,-1.1164e+04_r8, 9.8953e+04_r8,      &
       6.3735e+03_r8, 3.7240e+03_r8, 2.3588e+04_r8, 2.0106e+03_r8, 6.0596e+04_r8,      &
       4.6732e+04_r8,-3.1152e+05_r8,-1.7614e+04_r8, 9.4606e+04_r8, 3.4489e+03_r8,      &
       3.7961e+03_r8, 2.6755e+04_r8, 2.3096e+03_r8, 7.7763e+04_r8, 5.7709e+04_r8,      &
      -3.5374e+05_r8,-2.4401e+04_r8, 1.0057e+05_r8, 8.7498e+02_r8/

   data ((b720(ix,js),js=1,9),ix=1,5)/                                  &
       2.8378e+03_r8, 0.0000e+00_r8, 7.4867e+03_r8, 8.7260e+04_r8, 4.9146e+04_r8,      &
      -2.8306e+05_r8, 1.0229e+04_r8, 2.1829e+04_r8, 2.9344e+04_r8, 3.0405e+03_r8,      &
       0.0000e+00_r8, 1.0134e+04_r8, 9.3386e+04_r8, 4.9056e+04_r8,-2.9980e+05_r8,      &
       8.3963e+03_r8, 3.1439e+04_r8, 2.7751e+04_r8, 3.4447e+03_r8, 0.0000e+00_r8,      &
       1.4768e+04_r8, 9.9677e+04_r8, 4.7199e+04_r8,-3.2228e+05_r8, 4.6847e+03_r8,      &
       3.5809e+04_r8, 2.8838e+04_r8, 3.9763e+03_r8, 0.0000e+00_r8, 2.2241e+04_r8,      &
       1.0162e+05_r8, 3.8195e+04_r8,-3.2951e+05_r8,-3.5761e+02_r8, 3.1367e+04_r8,      &
       3.1066e+04_r8, 4.4566e+03_r8, 0.0000e+00_r8, 3.1280e+04_r8, 9.7135e+04_r8,      &
       2.7683e+04_r8,-3.2113e+05_r8,-4.8036e+03_r8, 2.4666e+04_r8, 3.2099e+04_r8/

   data ((b720(ix,js),js=1,9),ix=6,10)/                                 &
       1.9920e+03_r8, 4.3042e+03_r8, 3.4960e+04_r8, 9.2325e+04_r8, 1.9550e+04_r8,      &
      -3.0863e+05_r8,-7.8736e+03_r8, 2.0414e+04_r8, 3.1867e+04_r8, 2.3607e+03_r8,      &
       6.0047e+03_r8, 3.4158e+04_r8, 9.0097e+04_r8, 1.6570e+04_r8,-3.0384e+05_r8,      &
      -9.9396e+03_r8, 2.0701e+04_r8, 3.0825e+04_r8, 2.7001e+03_r8, 8.1576e+03_r8,      &
       3.0612e+04_r8, 8.9617e+04_r8, 1.4401e+04_r8,-2.9982e+05_r8,-1.1977e+04_r8,      &
       2.2595e+04_r8, 2.9273e+04_r8, 3.0360e+03_r8, 1.1081e+04_r8, 2.7697e+04_r8,      &
       9.1658e+04_r8, 1.1561e+04_r8,-3.0089e+05_r8,-1.4403e+04_r8, 2.5524e+04_r8,      &
       2.8020e+04_r8, 3.4215e+03_r8, 1.4885e+04_r8, 2.6879e+04_r8, 9.6042e+04_r8,      &   
       7.9992e+03_r8,-3.0885e+05_r8,-1.7495e+04_r8, 2.8632e+04_r8, 2.7565e+04_r8/

   data ((b720(ix,js),js=1,9),ix=11,15)/                                &
       3.8631e+03_r8, 1.9005e+04_r8, 2.7888e+04_r8, 1.0225e+05_r8, 4.8521e+03_r8,      &
      -3.2437e+05_r8,-2.2097e+04_r8, 3.3740e+04_r8, 2.7150e+04_r8, 4.3522e+03_r8,      &  
       2.2605e+04_r8, 3.0259e+04_r8, 1.0826e+05_r8, 1.4579e+03_r8,-3.4145e+05_r8,      &
      -2.7813e+04_r8, 3.9189e+04_r8, 2.6718e+04_r8, 4.8734e+03_r8, 2.5856e+04_r8,      &
       3.3158e+04_r8, 1.1340e+05_r8,-2.1817e+03_r8,-3.5946e+05_r8,-3.3255e+04_r8,      &
       4.5530e+04_r8, 2.5944e+04_r8, 5.3812e+03_r8, 2.8555e+04_r8, 3.5608e+04_r8,      &
       1.1671e+05_r8,-5.6934e+03_r8,-3.7309e+05_r8,-3.8155e+04_r8, 5.1594e+04_r8,      &
       2.4865e+04_r8, 5.8581e+03_r8, 3.1107e+04_r8, 3.6904e+04_r8, 1.1784e+05_r8,      &
      -9.1119e+03_r8,-3.7911e+05_r8,-4.2332e+04_r8, 5.6614e+04_r8, 2.3383e+04_r8/

   data ((b720(ix,js),js=1,9),ix=16,20)/                                &
       6.3262e+03_r8, 3.3720e+04_r8, 3.7274e+04_r8, 1.1891e+05_r8,-1.1690e+04_r8,      &
      -3.8488e+05_r8,-4.5891e+04_r8, 6.1655e+04_r8, 2.1870e+04_r8, 6.8035e+03_r8,      &
       3.6904e+04_r8, 3.7600e+04_r8, 1.2174e+05_r8,-1.2075e+04_r8,-3.9908e+05_r8,      &
      -4.8826e+04_r8, 6.8095e+04_r8, 2.1052e+04_r8, 7.2152e+03_r8, 4.0328e+04_r8,      &  
       3.8374e+04_r8, 1.2712e+05_r8,-1.0373e+04_r8,-4.2320e+05_r8,-5.1740e+04_r8,      &
       7.6305e+04_r8, 2.0951e+04_r8, 9.8640e+03_r8, 3.7426e+04_r8, 4.2372e+04_r8,      &
       1.3471e+05_r8,-6.8135e+03_r8,-4.5532e+05_r8,-5.4983e+04_r8, 8.5856e+04_r8,      &
       2.1689e+04_r8, 1.3560e+04_r8, 3.3322e+04_r8, 4.7176e+04_r8, 1.4419e+05_r8,      &
      -1.7431e+03_r8,-4.9373e+05_r8,-5.8662e+04_r8, 9.5890e+04_r8, 2.3540e+04_r8/

   data ((b720(ix,js),js=1,9),ix=21,25)/                                &
       1.8962e+04_r8, 2.8655e+04_r8, 5.2627e+04_r8, 1.5648e+05_r8, 4.6977e+03_r8,      &
      -5.4094e+05_r8,-6.3178e+04_r8, 1.0669e+05_r8, 2.6814e+04_r8, 2.6434e+04_r8,      &
       2.4076e+04_r8, 5.8634e+04_r8, 1.7285e+05_r8, 1.2578e+04_r8,-5.9982e+05_r8,      &
      -6.8826e+04_r8, 1.1902e+05_r8, 3.1514e+04_r8, 3.5147e+04_r8, 2.0130e+04_r8,      &
       6.5020e+04_r8, 1.9328e+05_r8, 2.1789e+04_r8,-6.6785e+05_r8,-7.5631e+04_r8,      &
       1.3244e+05_r8, 3.7235e+04_r8, 4.1827e+04_r8, 1.7320e+04_r8, 7.0904e+04_r8,      &
       2.1618e+05_r8, 3.0623e+04_r8,-7.3675e+05_r8,-8.2378e+04_r8, 1.4290e+05_r8,      &
       4.5585e+04_r8, 4.6294e+04_r8, 1.6071e+04_r8, 7.8131e+04_r8, 2.4682e+05_r8,      &  
       3.8078e+04_r8,-8.2274e+05_r8,-9.1157e+04_r8, 1.5361e+05_r8, 5.8726e+04_r8/

   data ((b720(ix,js),js=1,9),ix=26,31)/                                &
       4.9878e+04_r8, 1.6065e+04_r8, 8.8228e+04_r8, 2.8824e+05_r8, 4.2519e+04_r8,      &
      -9.3298e+05_r8,-1.0344e+05_r8, 1.6515e+05_r8, 7.8533e+04_r8, 5.3229e+04_r8,      &
       1.7740e+04_r8, 1.0255e+05_r8, 3.4331e+05_r8, 4.1675e+04_r8,-1.0732e+06_r8,      &
      -1.2083e+05_r8, 1.7696e+05_r8, 1.0553e+05_r8, 5.6585e+04_r8, 2.1720e+04_r8,      & 
       1.2220e+05_r8, 4.1397e+05_r8, 3.3029e+04_r8,-1.2465e+06_r8,-1.4427e+05_r8,      &
       1.8802e+05_r8, 1.3975e+05_r8, 5.9316e+04_r8, 2.8288e+04_r8, 1.4733e+05_r8,      &
       4.8733e+05_r8, 1.5630e+04_r8,-1.4301e+06_r8,-1.7016e+05_r8, 1.9433e+05_r8,      &
       1.7788e+05_r8, 6.2252e+04_r8, 3.8013e+04_r8, 1.8011e+05_r8, 5.5712e+05_r8,      &
      -1.8402e+04_r8,-1.6194e+06_r8,-1.9622e+05_r8, 2.0377e+05_r8, 2.1356e+05_r8,      &
       6.5578e+04_r8, 5.1423e+04_r8, 2.1882e+05_r8, 5.6668e+05_r8,-5.8779e+04_r8,      &
      -1.7453e+06_r8,-2.1194e+05_r8, 2.1181e+05_r8, 2.3210e+05_r8/

   data ((b720(ix,js),js=1,9),ix=32,37)/                                &
       6.9606e+04_r8, 6.8196e+04_r8, 2.6312e+05_r8, 5.9063e+05_r8,-1.5744e+05_r8,      &
      -1.8524e+06_r8,-2.3319e+05_r8, 2.6085e+05_r8, 2.1170e+05_r8, 7.5218e+04_r8,      &
       9.4375e+04_r8, 3.1888e+05_r8, 5.8441e+05_r8,-1.8982e+05_r8,-2.0276e+06_r8,      &
      -2.8095e+05_r8, 3.2730e+05_r8, 1.5852e+05_r8, 8.1818e+04_r8, 1.2971e+05_r8,      &
       3.6827e+05_r8, 6.1914e+05_r8,-2.8165e+05_r8,-2.2420e+06_r8,-3.4355e+05_r8,      &
       4.3737e+05_r8, 1.0553e+05_r8, 8.9628e+04_r8, 1.7800e+05_r8, 4.0374e+05_r8,      &
       5.5218e+05_r8,-2.6089e+05_r8,-2.4102e+06_r8,-4.0173e+05_r8, 4.6432e+05_r8,      &
       8.6117e+04_r8, 1.0036e+05_r8, 2.4647e+05_r8, 3.7601e+05_r8, 4.9780e+05_r8,      &
      -2.7995e+05_r8,-2.6053e+06_r8,-4.2444e+05_r8, 5.1942e+05_r8, 1.1274e+05_r8,      &
       1.1449e+05_r8, 3.4157e+05_r8, 3.5045e+05_r8, 4.8110e+05_r8,-2.9572e+05_r8,      &
      -2.9462e+06_r8,-4.6715e+05_r8, 6.1840e+05_r8, 1.4715e+05_r8/

   data ((b720(ix,js),js=1,9),ix=38,43)/                                &
       1.3273e+05_r8, 4.7433e+05_r8, 3.6299e+05_r8, 5.2029e+05_r8,-3.8660e+05_r8,      &    
      -3.6364e+06_r8,-5.3004e+05_r8, 9.1163e+05_r8, 1.8527e+05_r8, 1.5094e+05_r8,      &
       6.1556e+05_r8, 3.4710e+05_r8, 5.6221e+05_r8,-4.2073e+05_r8,-4.2353e+06_r8,      &
      -5.8904e+05_r8, 1.1559e+06_r8, 2.0012e+05_r8, 1.6568e+05_r8, 7.1379e+05_r8,      &
       2.4118e+05_r8, 6.1735e+05_r8,-3.9948e+05_r8,-4.6347e+06_r8,-5.9224e+05_r8,      &
       1.3508e+06_r8, 1.8977e+05_r8, 1.6599e+05_r8, 7.6711e+05_r8, 1.5253e+05_r8,      &
       7.5990e+05_r8,-2.0122e+05_r8,-4.9316e+06_r8,-6.2692e+05_r8, 1.3721e+06_r8,      &
       1.2698e+05_r8, 1.6400e+05_r8, 7.6327e+05_r8, 1.1905e+05_r8, 9.5620e+05_r8,      &
      -5.0089e+04_r8,-5.0316e+06_r8,-7.0061e+05_r8, 1.2124e+06_r8, 6.6342e+04_r8,      &
       1.6539e+05_r8, 8.6958e+05_r8, 9.0836e+04_r8, 1.1753e+06_r8, 8.8077e+04_r8,      &
      -5.4751e+06_r8,-7.6386e+05_r8, 1.2589e+06_r8, 2.1098e+04_r8/

      real(r8) :: o3pxfac(pver)      ! o3p cooling masking factors on WACCM vertical grids

      logical :: apply_co2_limit = .false.
      integer :: k1mb = 0
      
!================================================================================================
contains
!================================================================================================

  subroutine nlte_fomichev_init ( co2_mwi, n2_mwi, o1_mwi, o2_mwi, o3_mwi, no_mwi, apply_co2_limit_in )
     use interpolate_data, only : lininterp
     use ref_pres,         only : pref_edge, pref_mid

!     
!     Original version from Ray Roble
!     First adapted to CCM by F. Sassi - November 1999
!     
!---------------------------------------------------------------------------------
!     input:
!     RCO2 - CO2 volume mixing in the region below x=12.5
!     initial data from BLOCK DATA PCO2O3 come through common blocks
      
!     output: parameterization coefficients for both, the matrix parameterization 
!     (is used between x=2 and 12.5) and reccurence formula. 
!     AMAT,BMAT(43,9) - coefficients for the matrix parameterization

!-----------------------------------------------------------------

    real(r8), intent(in) :: o1_mwi                      ! O molecular weight
    real(r8), intent(in) :: o2_mwi                      ! O2 molecular weight
    real(r8), intent(in) :: o3_mwi                      ! O3 molecular weight
    real(r8), intent(in) :: co2_mwi                     ! CO2 molecular weight
    real(r8), intent(in) :: n2_mwi                      ! N2 molecular weight
    real(r8), intent(in) :: no_mwi                      ! NO molecular weight
    logical,  intent(in) :: apply_co2_limit_in

    real(r8), parameter :: p0=5.e-5_r8 ! TIE-GCM reference pressure in Pa            
    integer, parameter :: tgcmlevs = 29
    real(r8) :: pz(tgcmlevs)           ! TIE-GCM pressure grids (single resolution,pz=-7...7,dpz=0.5 ), dimensionless 
    real(r8) :: pp(tgcmlevs)           ! convert pz to Pascal  
    real(r8) :: xfac0(tgcmlevs)        ! masking factors on TIE-GCM pressure grids
    integer :: k

    apply_co2_limit = apply_co2_limit_in
    find_k1mb: do k = 1, pverp
       ! Find 1 mbar (or 100 Pa) level.
       if (pref_edge(k) > 100._r8) then
          k1mb = k
          exit find_k1mb
       endif
    end do find_k1mb
    if (masterproc) then
       write(iulog,'(a,l8)')  'nlte_fomichev_init: apply_co2_limit: ',apply_co2_limit
       write(iulog,'(a,i6,g12.6)') 'nlte_fomichev_init: check CO2 mixing ratios above 1-mbar level: ',k1mb,pref_mid(k1mb)
    endif
    
! set molecular weights
    co2_mw = co2_mwi
    n2_mw = n2_mwi
    o1_mw = o1_mwi
    o2_mw = o2_mwi
    o3_mw = o3_mwi
    no_mw = no_mwi

    ktop_co2cool = 1
    do k=1,pver
       if (pref_mid(k) < ptop_co2cool) ktop_co2cool = k
    enddo

    ! op3cooling masking factor  (from Kockarts and Peetermans [1981]
    xfac0(1:3)=.01_r8
    xfac0(4:10)=(/.05_r8,.1_r8,.2_r8,.4_r8,.55_r8,.7_r8,.75_r8/)
    xfac0(11:tgcmlevs) = .8_r8

    ! convert TIE-GCM pressure grid to Pascal

    pz(1)=-7.0_r8
    do k=2,tgcmlevs
       pz(k)=pz(k-1)+0.5_r8
    enddo
    do k=1,tgcmlevs
       pp(k)=p0*exp(-pz(k))
    enddo
    call lininterp( xfac0(tgcmlevs:1:-1), pp(tgcmlevs:1:-1), tgcmlevs, o3pxfac, pref_mid, pver )
    do k=1,pver
       if (pref_mid(k) > pp(1)) then
          o3pxfac(k)=0._r8
       else if (pref_mid(k) <= pp(tgcmlevs)) then
          o3pxfac(k)=xfac0(tgcmlevs)
       endif
    enddo

  end subroutine nlte_fomichev_init

  subroutine set_matrices( amat, bmat )

    implicit none
! calculate coefficients for the matrix paramerization:

    real(r8), intent(out) :: amat(nrfmlte,nrfmltelv)
    real(r8), intent(out) :: bmat(nrfmlte,nrfmltelv)

!-----------------------------------------------------------------
! Local vars:
    real(r8) :: rco2
    real(r8) ::  co2int(4), a
    integer :: i,j,isgn

    rco2 = chem_surfvals_get('CO2VMR')

    amat(1:nrfmlte,1:nrfmltelv)=0.0_r8
    bmat(1:nrfmlte,1:nrfmltelv)=0.0_r8
    do i = 1,nrfmlte
       do j = 1,nrfmltelv

          if((i.le.5).and.(j.eq.2)) goto 1
          isgn = int(sign(1._r8,a150(i,j))+sign(1._r8,a360(i,j))+          &
               sign(1._r8,a540(i,j))+sign(1._r8,a720(i,j)))
          co2int(1)=a150(i,j)/co2o(1)
          co2int(2)=a360(i,j)/co2o(2)
          co2int(3)=a540(i,j)/co2o(3)
          co2int(4)=a720(i,j)/co2o(4)
          if(isgn.eq.-4) then
             co2int(1) = log(-co2int(1))
             co2int(2) = log(-co2int(2))
             co2int(3) = log(-co2int(3))
             co2int(4) = log(-co2int(4))
             a = -exp(a18lin(rco2,co2o,co2int,1,4))
          else if (isgn.eq.4) then
             co2int(1) = log(co2int(1))
             co2int(2) = log(co2int(2))
             co2int(3) = log(co2int(3))
             co2int(4) = log(co2int(4))
             a = exp(a18lin(rco2,co2o,co2int,1,4))
          else
             call a18int(co2o,co2int,rco2,a,4,1)
          end if
          amat(i,j)=a*rco2
          
          isgn = int(sign(1._r8,b150(i,j))+sign(1._r8,b360(i,j))+          &
               sign(1._r8,b540(i,j))+sign(1._r8,b720(i,j)))
          co2int(1)=b150(i,j)/co2o(1)
          co2int(2)=b360(i,j)/co2o(2)
          co2int(3)=b540(i,j)/co2o(3)
          co2int(4)=b720(i,j)/co2o(4)
          if(isgn.eq.-4) then
             co2int(1) = log(-co2int(1))
             co2int(2) = log(-co2int(2))
             co2int(3) = log(-co2int(3))
             co2int(4) = log(-co2int(4))
             a = -exp(a18lin(rco2,co2o,co2int,1,4))
          else if (isgn.eq.4) then
             co2int(1) = log(co2int(1))
             co2int(2) = log(co2int(2))
             co2int(3) = log(co2int(3))
             co2int(4) = log(co2int(4))
             a = exp(a18lin(rco2,co2o,co2int,1,4))
          else
             call a18int(co2o,co2int,rco2,a,4,1)
          end if
          bmat(i,j)=a*rco2
1         continue
       enddo
    enddo

    return
  endsubroutine set_matrices


!==================================================================================================

      subroutine nlte_fomichev_calc (lchnk,ncol,pmid,pint,t,xo2,xo,xo3,xn2,xco2,coolf,&
                                     co2cool_out, o3cool_out, c2scool_out )
         use time_manager,  only: get_nstep

!
!     author: F. Sassi (Dec, 1999)
!
!------------------------------------------------------------------
!
!     This is  routine prep arrays to be passed  to Fomichev' scheme
!
!     RADFMCINTI should have been called at the beginning of the run to
!     create arrays used in this parameterization.
!
!     Concentrations in input are expected in mass mixing ratios 
!     and are converted to volume mixing ratios
!
!     Because Fomichev scheme has its own grid which runs from 
!     ground upward, arrays are prepared with vertical indexing
!     inverted with respect to the CCM convention. However,
!     arrays in input are expected in the standard convention
!     (i.e., top is level 1).
!
!     The pressures at mid-points and interfaces need to be known 
!     Conversion to normalized X coordinate is carried out here.
!
!     Cooling rates are calculated in units of 
!
!                      dT
!           COOLF = Cp ---
!                      dt
!
!     where Cp is the specific heat (ergs g^-1 K^-1), dT/dt is the 
!     temperature tendency (K s^-1). Therefore, units of cooling
!     are 
!     
!         (erg g^-1 K^-1) (K s^-1) = cm^2 s^-3
!
!     COOLF is converted to J/kg/s before output. 
!
!
!*****************************************************************
!
!-----------------------------------------------------------------
implicit none
!-----------------------------------------------------------------

!     Input variables
      integer, intent(in) :: ncol                          ! number of atmospheric columns
      integer, intent(in) :: lchnk                         ! chunk identifier

      real(r8), intent(in) :: pmid(pcols,pver)             ! model pressure at mid-point
      real(r8), intent(in) :: pint(pcols,pverp)            ! model pressure at interfaces
      real(r8), intent(in) :: t(pcols,pver)                ! Neutral temperature (K)
      real(r8), intent(in) :: xco2(pcols,pver)             ! CO2 profile
      real(r8), intent(in) :: xn2(pcols,pver)              ! N2 profile
      real(r8), intent(in) :: xo3(pcols,pver)              ! O3 profile
      real(r8), intent(in) :: xo(pcols,pver)               ! O profile
      real(r8), intent(in) :: xo2(pcols,pver)              ! O2 profile

!     Output variables
      real(r8), intent(out) :: coolf(pcols,pver)          ! Total cooling
      real(r8), intent(out) :: co2cool_out(pcols,pver)        ! CO2 cooling 
      real(r8), intent(out) :: o3cool_out(pcols,pver)         ! O3 cooling
      real(r8), intent(out) :: c2scool_out(pcols,pver)        ! Cooling to Space

!    Local variables
      
      real(r8) rmo2                         ! O2 molecular weight
      real(r8) rmo                          ! O molecular weight
      real(r8) rmn2                         ! N2 molecular weight
      real(r8) rmco2                        ! CO2 molecular weight
      real(r8) rmo3                         ! O3 molecular weight
      real(r8) xnorm(pcols,pver)            ! normalized X p.s.h. at midpoints
      real(r8) xnori(pcols,pverp)           ! normalized X p.s.h. at interfaces
      real(r8) dxnorm(pcols,pver)           ! xnorm(k+1)-xnorm(k)
      real(r8) presm(pcols,pver)            ! pressure at midpoint (dyn/cm^2)
      real(r8) presi(pcols,pverp)           ! pressure at interfaces (dyn/cm^2)
      real(r8) dpi(pcols,pver)              ! pressure diff. between interfaces (dyn/cm^2)
      real(r8) mwair(pcols,pver)            ! mean air molecular weight (g/mole)
      real(r8) ndenair(pcols,pver)          ! mean air number density (cm**-3)
      real(r8) colco2(pcols,pver)           ! CO2 column number density
      real(r8) uco2(pcols,nrfm)             ! column CO2
      real(r8) dummyg(pver)                 ! dummy
      real(r8) dummyx(pver)                 ! dummy
      real(r8) dummyf(nrfm)                 ! dummy
      real(r8) hco2(pcols,nrfm)             ! CO2 cooling in Fomichev grid
      real(r8) ho3(pcols,nrfm)              ! O3 cooling in Fomichev grid
      real(r8) tf(pcols,nrfm)               ! neutral temp interpolated to Fomichev grid
      real(r8) vn2f(pcols,nrfm)             ! N2 vmr interpolated to Fomichev grid
      real(r8) vo3f(pcols,nrfm)             ! O3 vmr interpolated to Fomichev grid
      real(r8) vof(pcols,nrfm)              ! O vmr interpolated to Fomichev grid
      real(r8) vco2f(pcols,nrfm)            ! CO2 vmr interpolated to Fomichev grid
      real(r8) vo2f(pcols,nrfm)             ! O2 vmr interpolated to Fomichev grid
      real(r8) mwairf(pcols,nrfm)           ! Mean air molecular weight interpolated to Fomichev grid
      real(r8) ndenf(pcols,nrfm)            ! Mean air no. density interpolated to Fomichev grid
      real(r8) flux(pcols)                  ! Flux boundary condition for cool-to-space
      real(r8) vo2(pcols,pver)              ! O2 vmr
      real(r8) vo(pcols,pver)               ! O vmr
      real(r8) vo3(pcols,pver)              ! O3 vmr
      real(r8) vn2(pcols,pver)              ! N2 vmr
      real(r8) vco2(pcols,pver)             ! CO2 vmr
      real(r8) co2cooln(pcols,pver)         ! CO2 cooling 
      real(r8) o3cooln(pcols,pver)          ! O3 cooling
      real(r8) hc2s(pcols,pver)             ! cool to space heating
      real(r8) ps0                          ! Reference (surface) pressure
      real(r8) ti(pcols,pver)               ! T(PVER:1:-1)
      real(r8) alam(pcols,nrfm)             ! LAMBDA
      real(r8) djm(pcols,nrfm)              ! DJM in recurrence formula
      real(r8) dj0(pcols,nrfm)              ! DJ0 in recurrence formula
      real(r8) aajm(pcols,nrfm)             ! AAJM in recurrrence formula
      real(r8) aaj0(pcols,nrfm)             ! AJJ0 in recurrence formula
      real(r8) zhgt(pcols,pver)             ! approx. elevation in cm
      real(r8) grav(pcols,pver)             ! accelration of gravity in cm/s^2
      real(r8) :: wrk(pcols)
      
      integer k
      integer i
      integer kinv                     ! inverted vertical index (bottom up)
      real(r8), parameter :: co2_limit = 720.e-6_r8
      integer :: nstep
      character(len=200) :: errmsg
      real(r8) :: latdeg, londeg

!----------------------------------------------------------------

!     Define molecular weights
      rmco2 = co2_mw
      rmo3  = o3_mw
      rmo2  = o2_mw
      rmo   = o1_mw
      rmn2  = n2_mw


      coolf(1:ncol,1:pver)=0.0_r8

      arad=rearth*1e2_r8  ! initialize planet's radius (cm)

!-----------------------------------------------------------------
!     The pressure at mid point is converted to normalized x coordinate
!                  X = LN (P0 / PRES)
!     where P0 = 1e6 dyn/cm^2 ; PRES is in dyn/cm^2
!     Note: to convert from Pa to dyn/cm^2, multyply by 10
!-----------------------------------------------------------------
      ps0=1.000e6_r8
      do k=1,pver
         kinv = pver-k+1
         do i=1,ncol
            presm(i,k) = pmid(i,kinv)*10._r8
            xnorm(i,k) = log(ps0/presm(i,k))
!     Calculate pressure at interfaces
            presi(i,k) = pint(i,kinv+1)*10._r8
            xnori(i,k) = log(ps0/presi(i,k))
         enddo
      enddo
      presi(:ncol,pverp) = pint(:ncol,1)*10._r8
      xnori(:ncol,pverp) = log(ps0/presi(:ncol,pverp))

!-----------------------------------------------------------------
!     Calculate layer thikcness (DPI).
!     For each pressure interface the following is true:
!     
!            Pint(k) = P0 * exp (-Xint(k))
!
!     Thus,
!
!            DPI(k)= Pint(k)-Pint(k+1)= P0 * exp (-Xint(k)) * DX
!
!     where,
!
!            DX = 1 - exp ( - (Xint(k+1)-Xint(k)) )
!
!-----------------------------------------------------------------
      do k=pver,1,-1
         do i=1,ncol
            dxnorm(i,k)=xnori(i,k+1)-xnori(i,k)
         enddo
      enddo

!     Pressure difference between interfaces (positive downward)
      do k=1,pver
         do i=1,ncol
            dpi(i,k)=ps0*exp(-xnorm(i,k))*(1._r8-exp(-dxnorm(i,k)))
         enddo
      enddo

!-----------------------------------------------------------------
!     Calculate molecular weight (g/mol) of mean air MWAIR:
!      
!               MWAIR= 1. / [ Sum(i) MMR(i)/MW(i) ] 
!
!     where MMR(i) are  mass mixing ratio of O2,
!     O, N2, and MW(i) are the corresponding
!     molecular weights.
!-----------------------------------------------------------------
      mwair(1:ncol,1:pver)=0.0_r8
      do k=1,pver
         kinv=pver-k+1
         do i=1,ncol
            mwair(i,k)=1._r8/ (                         &
                 xo2(i,kinv)/rmo2                    &
                +xo(i,kinv)/rmo                      &
                +xn2(i,kinv)/rmn2                    &
                +xo3(i,kinv)/rmo3                    &
                +xco2(i,kinv)/rmco2                  &
                )
         enddo
      enddo

!
!     Elevation is calculated via integration of 
!     the hydrostatic equation 
!
!
!           Sum(k) g(k) dz = Sum(k) PHI(k)
!
!     where 
!                     g0*a^2
!               g(k)= ------
!                     (a+z)^2
!
!     where g0=980 cm/s^2;a=6.37e8
!
!     and
!
!                                UR * T(i)  DPI(i)
!            PHI(k) = Sum(i=1,k)  ------    ----
!                                  MW(i)    P(i)
!     
!     where UR is the gas universal constant, T is temperature,
!     MW is the mean air molecular weight, DPI is the pressure 
!     layer thickness and P is the mid-point pressure
!     Then,
!
!                   PHI * a
!              z = ---------
!                  (a*g0-PHI)
!
!     do i=1,ncol
!        phi=0.0
!        do k=1,pver
!           kinv=pver-k+1
!           phi=phi+ur*t(i,kinv)/(mwair(i,k))*dpi(i,k)/presm(i,k)
!           zhgt(i,k)=phi*arad/(arad*grav0-phi)
!           grav(i,k)=grav0*(arad/(arad+zhgt(i,k)))**2
!        enddo
!     enddo

      wrk(:ncol) = 0._r8

      do k=1,pver
         kinv=pver-k+1
         do i=1,ncol
            wrk(i) = wrk(i) + ur*t(i,kinv)/(mwair(i,k))*dpi(i,k)/presm(i,k)
            zhgt(i,k) = wrk(i)*arad/(arad*grav0 - wrk(i))
            grav(i,k) = grav0*(arad/(arad + zhgt(i,k)))**2
         enddo
      enddo

      do k = 1,pver
         kinv=pver-k+1
         do i=1,ncol
!-----------------------------------------------------------------
!     Convert mmr to vmr
!-----------------------------------------------------------------
            vo2 (i,k) = xo2 (i,kinv) *mwair(i,k)/rmo2
            vo  (i,k) = xo  (i,kinv) *mwair(i,k)/rmo
            vo3 (i,k) = xo3 (i,kinv) *mwair(i,k)/rmo3
            vn2 (i,k) = xn2 (i,kinv) *mwair(i,k)/rmn2
            vco2(i,k) = xco2(i,kinv) *mwair(i,k)/rmco2

! CGB - The Formichev scheme was not designed to support CO2 > 720 ppmv, so
! limit the amount of CO2 used to 720 ppmv. This keeps the model stable, but
! may yield an incorrect scientific result. It would be nice to extend this
! routine to support higher CO2 values. Putting the limiter here means that
! that the other constituents will have their proper mixing ratio caclulated
! (i.e. the mwair is correct), but vco2 will be limited.  Abort the run if CO2
! exceeds the limit at altitudes above 1 mbar unless apply_co2_limit=.true.

            if ((vco2(i,k)>co2_limit).and.(.not.apply_co2_limit).and.(kinv<k1mb)) then
               nstep = get_nstep()
               latdeg = get_rlat_p(lchnk,i)*180._r8/pi
               londeg = get_rlon_p(lchnk,i)*180._r8/pi
               write(errmsg,fmt='(a,i12,2(i6),g12.4,2(f8.2),g12.4)') &
                    'nlte_fomichev_calc: CO2 has exceeded the limit: nstep,i,k,press(Pa),lon,lat,vco2(vmr)=',&
                    nstep,i,kinv, pmid(i,kinv), londeg, latdeg, vco2(i,k)
               write(iulog,*) trim(errmsg)
               call endrun(trim(errmsg))
            endif

!-----------------------------------------------------------------
!     Calculate mean air number density ( cm^(-3) )
!
!                              P(k)
!                    n  = -------------
!                              Kb*T(k)
!     where:
!           P(k)    is pressure at mid-point
!           AKBL    is Boltzman constant
!           T       is neutral temperature
!-----------------------------------------------------------------
            ndenair(i,k) = presm(i,k)/(akbl*t(i,kinv))
         enddo

      enddo

      ! apply the CO2 limiter for all levels and columns
      where ( vco2(:ncol,:) > co2_limit )
         vco2(:ncol,:) = co2_limit
      end where
      
!-----------------------------------------------------------------
!     Calculate CO2 vertical column above each level
!     
!     At each mid-point vertical level (j), the following sum is calculated
!
!                                         
!            COLCO2(j) = Sum_(i=ztop:j:-1) NDENAIR(I)*VCO2(I)*DZ = ...
!                      
!                                            ANAV * VCO2(I)
!                      = Sum_(i=pver:j:-1)  ---------------- * DP
!                                            GRAV * MWAIR(I)
!
!     where ANAV is the Avogadro no., VCO2 is CO2 vmr, GRAV is the
!     accelaration of gravity, MWAIR is mean molecular weight, DP
!     is the pressure increment downward.
!     As boundary condition at PVER, it is assumed that VCO2 
!     and NDENAIR stay constant above that level.
!-----------------------------------------------------------------
!     do i=1,ncol
!        colco2(i,pver)=anav/grav(i,pver)*vco2(i,pver)/mwair(i,pver)*dpi(i,pver)
!        do k=pver-1,1,-1
!           colco2(i,k)=colco2(i,k+1)                                       &
!               +anav/grav(i,k)*vco2(i,k)/mwair(i,k)*dpi(i,k)
!        enddo
!     enddo

      colco2(:ncol,pver) = anav/grav(:ncol,pver)*vco2(:ncol,pver) &
			   /mwair(:ncol,pver)*dpi(:ncol,pver)
      do k=pver-1,1,-1
         do i=1,ncol
            colco2(i,k)=colco2(i,k+1)                                       &
                +anav/grav(i,k)*vco2(i,k)/mwair(i,k)*dpi(i,k)
         enddo
      enddo

!-----------------------------------------------------------------
!     Linear interpolation from CCM grid defined in XNORM(1:PVER)
!     to Fomichev grid defined in XR(1:NRFM)
!     Interpolation is carried out using mod PRFLINV which
!     is adapted from A18LINV (originally written for TIME/GCM).
!     All fields are interpolated. If XR levels are beyond the
!     limit of the actual CCM grid, zero values are inserted in
!     the interpolated arrays. Proper handling of the arrays is done
!     in VICCOOLN.
!-----------------------------------------------------------------

!-----------------------------------------------------------------
!     Create TI array which contains T(PVER:1:-1)
!-----------------------------------------------------------------
      ti(1:ncol,1:pver)=t(1:ncol,pver:1:-1)

      do i=1,ncol

         dummyx(1:pver)=xnorm(i,1:pver)

!     Temperature
         dummyg(1:pver)=ti(i,1:pver)
         call a18linvne (xr,dummyf,dummyx,dummyg,pver,nrfm)
         tf(i,1:nrfm)=dummyf(1:nrfm)

!     O3
         dummyg(1:pver)=vo3(i,1:pver)
         call a18linvne (xr,dummyf,dummyx,dummyg,pver,nrfm)
         vo3f(i,1:nrfm)=dummyf(1:nrfm)

!     O2
         dummyg(1:pver)=vo2(i,1:pver)
         call a18linvne (xr,dummyf,dummyx,dummyg,pver,nrfm)
         vo2f(i,1:nrfm)=dummyf(1:nrfm)

!     N2
         dummyg(1:pver)=vn2(i,1:pver)
         call a18linvne (xr,dummyf,dummyx,dummyg,pver,nrfm)
         vn2f(i,1:nrfm)=dummyf(1:nrfm)

!     O
         dummyg(1:pver)=vo(i,1:pver)
         call a18linvne (xr,dummyf,dummyx,dummyg,pver,nrfm)
         vof(i,1:nrfm)=dummyf(1:nrfm)
         
!     CO2
         dummyg(1:pver)=vco2(i,1:pver)
         call a18linvne (xr,dummyf,dummyx,dummyg,pver,nrfm)
         vco2f(i,1:nrfm)=dummyf(1:nrfm)

!     COLCO2
         dummyg(1:pver)=colco2(i,1:pver)
         call a18linvne (xr,dummyf,dummyx,dummyg,pver,nrfm)
         uco2(i,1:nrfm)=dummyf(1:nrfm)

!     DEN
         dummyg(1:pver)=ndenair(i,1:pver)
         call a18linvne (xr,dummyf,dummyx,dummyg,pver,nrfm)
         ndenf(i,1:nrfm)=dummyf(1:nrfm)
         
!     AM
         dummyg(1:pver)=mwair(i,1:pver)
         call a18linvne (xr,dummyf,dummyx,dummyg,pver,nrfm)
         mwairf(i,1:nrfm)=dummyf(1:nrfm)
         
      enddo


!     Use recurrence relation to calculate AL coefficents
      call recur (ncol,uco2,tf,vn2f,vo2f,vof,ndenf,alam,djm,dj0,aajm,aaj0)


!     Do LTE and NLTE parts of cooling
      call viccooln (ncol,alam,djm,dj0,aajm,aaj0,tf,vco2f,vo3f,mwairf,flux,hco2,ho3)


!     Interpolate from Fomichev grid to CCM grid
      do i=1,ncol
         dummyx(1:pver)=xnorm(i,1:pver)

!     HCO2
         dummyf(1:nrfm)=hco2(i,1:nrfm)
         call a18linvne (dummyx,dummyg,xr,dummyf,nrfm,pver)
         co2cooln(i,1:pver)=dummyg(1:pver)

!     HO3
         dummyf(1:nrfm)=ho3(i,1:nrfm)
         call a18linvne (dummyx,dummyg,xr,dummyf,nrfm,pver)
         o3cooln(i,1:pver)=dummyg(1:pver)

      enddo

!     Do cool-to-space component of cooling
      call cool2space (ncol,ti,mwair,vn2,vo2,vo,vco2,ndenair,xnorm,flux,hc2s)

!     Calculate total cooling

      ! Above ptop_co2cool use cool to space approx.
      do k=1,ktop_co2cool-1
         kinv=pver-k+1
         do  i=1,ncol
            ! Convert to J/kg/s
            coolf(i,k) = (o3cooln(i,kinv) + hc2s(i,kinv)) * 1.e-4_r8
         enddo
      enddo
      ! Below ptop_co2cool use nlte calculation
      do k=ktop_co2cool,pver
         kinv=pver-k+1
         do  i=1,ncol
            ! Convert to J/kg/s
            coolf(i,k) = (co2cooln(i,kinv) + o3cooln(i,kinv)) * 1.e-4_r8
         enddo
      enddo

      ! diagnostics ...
      do k=1,pver
         kinv=pver-k+1
         co2cool_out(:ncol,k) = co2cooln(:ncol,kinv) * 1.e-4_r8
         o3cool_out(:ncol,k) = o3cooln(:ncol,kinv) * 1.e-4_r8
         c2scool_out(:ncol,k) = hc2s(:ncol,kinv) * 1.e-4_r8
      enddo

      return
    end subroutine nlte_fomichev_calc

!====================================================================================

  subroutine a18linvne (x,y,xn,yn,n,imax)                               

      implicit none

!     ****                                                               
!     ****     This procedure performs linear interpolation within the   
!     ****     table defined by the N points (XN(NN),Y(NN)).             
!     ****     Where:                                                    
!     ****                                                               
!     ****       NN = 1,N,1                                              
!     ****                                                               
!     ****       XN(NN) < XN(NN+1) for NN = 1,N-1            
!     ****                                                               
!     ****     Parameters:                                               
!     ****                                                               
!     ****       X(IMAX) = array of IMAX x-values at which linear        
!     ****                 interpolation is required                     
!     ****                                                               
!     ****       XN(N) = array of N abscissae at which function values   
!     ****               are given                                       
!     ****                                                               
!     ****       YN(N) = function values corresponding to abscissae,     
!     ****               XN(N)                                           
!     ****                                                               
!     ****     Output:                                                   
!     ****                                                               
!     ****      Y(IMAX)  The IMAX interpolated values are                
!     ****                     returned in this array                    
!     ****                                                               
!     
!     It has been modified as follows: 
!     if points X are outside the range X(1)..X(N), the last values
!     are assigned. That is
!     IF X(I) > XN(N) THEN Y(I)=YN(N)
!     IF X(I) < XN(1) THEN Y(I)=YN(1)
!     

!     Arguments
      integer, intent(in) :: imax
      integer, intent(in) :: n
      real(r8), intent(out) :: y(imax)
      real(r8), intent(in) :: x(imax)
      real(r8), intent(in) :: xn(n)
      real(r8), intent(in) :: yn(n)


!     Local variables
      integer kk(imax)                 
      integer nn
      integer i

!     ****                                                               
!     ****     Where:                                                    
!     ****       Y(IMAX) is vector output                                
!     ****                                                               
!     ****       KK is work space                                        
!     ****                                                               
!     ****                                                               
!     ****     Initialize array KK                                       
!     ****                                                               

      do i = 1,imax                                                      
         kk(i) = 0                                                        
      enddo                                                              

!     ****                                                               
!     ****     Locate interval in (XN,YN) in table containing X(I)       
!     ****                                                               

      do nn = 1,n-1                                                      
         do i = 1,imax                                                    
            kk(i) = merge(nn+1,kk(i),(xn(nn+1)-x(i))*(x(i)-xn(nn))>=0._r8)    
         enddo                                                            
      enddo                                                              

!     ****                                                               
!     ****     Check for                                                 
!     ****                                                               
!     ****       X(I) < XN(1),  X(I) > X(N)                              
!     ****                                                               
!     ****       and use linear extrapolation if necessary               
!     ****                                                               

      do i = 1,imax                                                      
         kk(i) = merge(-1,kk(i),xn(1)-x(i)>=0._r8)                            
         kk(i) = merge(-2,kk(i),x(i)-xn(n)>=0._r8)                            
      enddo                                                              

!     ****                                                               
!     ****     Perform interpolation prescribed above                    
!     ****                                                               

      y(:) = 0._r8
      
      do i = 1,imax      

         if (kk(i).gt.0) then

            y(i) = (                            &
                 yn(kk(i)-1)*(xn(kk(i))-x(i))   &
                 + yn(kk(i))*(x(i)-xn(kk(i)-1)) &
                 )/(xn(kk(i))-xn(kk(i)-1))         

         else if (kk(i).eq.-1) then

            y(i)=yn(1)

         else if (kk(i).eq.-2) then

            y(i)=yn(n)

         endif

      enddo                                                              

      return
 end subroutine  a18linvne


!=======================================================================================

      subroutine recur (ncol,uco2,tf,vn2f,vo2f,vof,ndenf,alam,djm,dj0,aajm,aaj0)

!     Originally written by R. Roble,  modified by F. Sassi (Nov., 1999)

!-----------------------------------------------------------------
implicit none
!-----------------------------------------------------------------


      integer, intent(in) :: ncol                          ! number of atmospheric columns
      
      real(r8), intent(in) :: UCO2(pcols,nrfm)
      real(r8), intent(in) :: tf(pcols,nrfm)
      real(r8), intent(in) :: vn2f(pcols,nrfm)
      real(r8), intent(in) :: vo2f(pcols,nrfm)
      real(r8), intent(in) :: vof(pcols,nrfm)
      real(r8), intent(in) :: ndenf(pcols,nrfm)

      real(r8), intent(out) :: alam(pcols,nrfm)
      real(r8), intent(out) :: djm(pcols,nrfm)
      real(r8), intent(out) :: dj0(pcols,nrfm)
      real(r8), intent(out) :: aajm(pcols,nrfm)
      real(r8), intent(out) :: aaj0(pcols,nrfm)

      real(r8) CO2INT(nrfmco2)
      real(r8) UREF(nrfmco2)
      real(r8) A(pcols)                  
      real(r8) COR(pcols)
      real(r8) UC(pcols)                                           
      real(r8) al(pcols,nrfm)

      real(r8) tt
      real(r8) y
      real(r8) zn2
      real(r8) zo2
      real(r8) zz
      real(r8) rko

      integer ks
      integer i
      integer k
      integer km


!****this constant should be moved to an intialization routine

!                                                                        
!                                                                        
!     ****  UCO2 (CO2 COLUMN AMOUNT) FOR CO2                         
!                                                                        
!                                                                        
!     **** CALCULATE COEFICIENTS FOR THE RECCURENCE FORMULA:             
!                                                                        
!     **** BETWEEN X=12.5 AND 13.75 THESE COEFFICIENTS (AL) ARE          
!     **** CALCULATED USING CORRECTIONS TO THE ESCAPE FUNCTION.          
!     **** STARTING FROM X=14.00 AND ABOVE THE PARAMETERIZATION          
!     **** COEFFICIENTS ARE EQUAL TO THE ESCAPE FUNCTION.                
!                        

      al(1:ncol,1:nrfm)=0.0_r8
      ks=0
      do k=1,nrfm                                                         
         
         if (xr(k).ge.12.5_r8 .and. xr(k).le.13.75_r8) then
            ks=ks+1

            co2int(1) = cor150(ks)
            co2int(2) = cor360(ks)
            co2int(3) = cor540(ks)
            co2int(4) = cor720(ks)
            uref(1) = uco2co(ks)*150._r8/360._r8
            uref(2) = uco2co(ks)
            uref(3) = uco2co(ks)*540._r8/360._r8
            uref(4) = uco2co(ks)*720._r8/360._r8
            do  i=1,ncol
               uc(i) = uco2(i,k)
            enddo
            call a18linv(uc,a,uco2ro,alo,51,ncol)
            call a18linv(uc,cor,uref,co2int,4,ncol)
            do i=1,ncol
               al(i,k) = exp(cor(i)+a(i))
            enddo

         endif
      enddo

      do k=1,nrfm

         if (xr(k).ge.14.00_r8) then

            do i=1,ncol
               uc(i) = uco2(i,k)
            enddo
            call a18linv(uc,a,uco2ro,alo,51,ncol)
            do  i=1,ncol
               al(i,k) = exp(a(i))
            enddo
            
         endif

      enddo

! 
!     Calculate ALAM
!                     
      alam(1:ncol,1:nrfm)=0.0_r8
      do  k=1,nrfm

!     ALAM is used only for p.s.h. >= 12.5
!     If the current level is below 12.5 s.h., then do nothing
         if (xr(k).ge.12.5_r8) then

            do i=1,ncol                                                  

!                                                                     
!     ****  CO2-O2 AND CO2-N2 V-T CONSTANTS                           
!                                                                     
               tt   = tf(i,k) 
               y    = tt**(-1._r8/3._r8)                                          
               zn2  = 5.5e-17_r8*sqrt(tt)+6.7e-10_r8*exp(-83.8_r8*y)            
               zo2  = 1.e-15_r8*exp(23.37_r8-230.9_r8*y+564._r8*y*y)            
               rko  = 3.0e-12_r8                                                    

!                                                                     
!     ****  COLLISIONAL DEACTIVATION RATE:                            
!                                                                     
               zz   = (vn2f(i,k)*zn2 +  vo2f(i,k)*zo2 +  vof (i,k)*rko)*ndenf(i,k)

!
!     ****
!
               alam(i,k) = a10/( a10+zz )

            enddo  ! end-loop in longitude
         
         endif

      enddo        ! end-loop in levels

!     Calculate coefficients of recurrence formula
!     This coefficients are used for 12.75=< p.s.h.=<16.5
!     Outside this range do nothing
!     It uses ALAM (p.s.h. >= 12.5)

      djm(1:ncol,1:nrfm)=0.0_r8
      dj0(1:ncol,1:nrfm)=0.0_r8
      aajm(1:ncol,1:nrfm)=0.0_r8
      aaj0(1:ncol,1:nrfm)=0.0_r8
      do k=1,nrfm

         if (xr(k).ge.12.75_r8 .and. xr(k).le.16.5_r8) then

            km=k-1

            do i=1,ncol

               djm(i,k)  = +.25_r8*(3._r8*al(i,km) +    al(i,k) )
               dj0(i,k)  = +.25_r8*(   al(i,km) + 3._r8*al(i,k) )

               aajm(i,k) = 1._r8-alam(i,km) * ( 1._r8-djm(i,k) )
               aaj0(i,k) = 1._r8-alam(i,k ) * ( 1._r8-dj0(i,k) )

            enddo    ! end-loop in longitude

         endif

      enddo          ! end-loop in levels

      return
  end subroutine recur
                                                                
!======================================================================================

      subroutine a18linv (x,y,xn,yn,n,imax)                               

      implicit none

!     ****                                                               
!     ****     This procedure performs linear interpolation within the   
!     ****     table defined by the N points (XN(NN),Y(NN)).             
!     ****     Where:                                                    
!     ****                                                               
!     ****       NN = 1,N,1                                              
!     ****                                                               
!     ****       XN(NN) < XN(NN+1) for NN = 1,N-1            
!     ****                                                               
!     ****     Parameters:                                               
!     ****                                                               
!     ****       X(IMAX) = array of IMAX x-values at which linear        
!     ****                 interpolation is required                     
!     ****                                                               
!     ****       XN(N) = array of N abscissae at which function values   
!     ****               are given                                       
!     ****                                                               
!     ****       YN(N) = function values corresponding to abscissae,     
!     ****               XN(N)                                           
!     ****                                                               
!     ****     Output:                                                   
!     ****                                                               
!     ****      Y(IMAX)  The IMAX interpolated values are                
!     ****                     returned in this array                    
!     ****                                                               

!     Input variables
      integer imax
      integer n
      real(r8) y(imax)
      real(r8) x(imax)
      real(r8) xn(n)
      real(r8) yn(n)

      integer i

!     Local variables
      integer KK(IMAX)                 
      integer NN

!     ****                                                               
!     ****     Where:                                                    
!     ****       Y(IMAX) is vector output                                
!     ****                                                               
!     ****       KK is work space                                        
!     ****                                                               
!     ****                                                               
!     ****     Initialize array KK                                       
!     ****                                                               

      do i = 1,imax                                                      
         kk(i) = 0                                                        
      enddo                                                              

!     ****                                                               
!     ****     Locate interval in (XN,YN) in table containing X(I)       
!     ****                                                               

      do nn = 1,n-1                                                      
         do i = 1,imax                                                    
            kk(i) = merge(nn+1,kk(i),(xn(nn+1)-x(i))*(x(i)-xn(nn))>=0._r8)    
         enddo                                                            
      enddo                                                              

!     ****                                                               
!     ****     Check for                                                 
!     ****                                                               
!     ****       X(I) < XN(1),  X(I) > X(N)                              
!     ****                                                               
!     ****       and use linear extrapolation if necessary               
!     ****                                                               

      do i = 1,imax                                                      
         kk(i) = merge(2,kk(i),xn(1)-x(i)>=0._r8)                            
         kk(i) = merge(n,kk(i),x(i)-xn(n)>=0._r8)                            
      enddo                                                              

!     ****                                                               
!     ****     Perform interpolation prescribed above                    
!     ****                                                               

      do i = 1,imax                                                      
         y(i) = (yn(kk(i)-1)*(xn(kk(i))-x(i)) + yn(kk(i))*    &
             (x(i)-xn(kk(i)-1)))/(xn(kk(i))-xn(kk(i)-1))         
      enddo                                                              

      return                                                             
 end  subroutine a18linv


!========================================================================================

      subroutine viccooln (ncol,alam,djm,dj0,aajm,aaj0,tv,co2,o3,am,flux,hco2,ho3 )

!
!     Original version from Ray Roble.
!     Adapted to CCM by F. Sassi (Nov. 1999)
!

!-----------------------------------------------------------------
implicit none
!-----------------------------------------------------------------

!
!     **** This is the mod that calculates LTE and NLTE components
!     **** of the cooling rates.
!

! XL(17)  - the parameters for NLTE region (12.5 <= X <= 16.5)    

!     Input variables

      integer, intent(in) :: ncol                          ! number of atmospheric columns

      real(r8), intent(in) :: tv(pcols,nrfm)               ! neutral temp interpolated to Fomichev grid
      real(r8), intent(in) :: o3(pcols,nrfm)               ! O3 vmr interpolated to Fomichev grid
      real(r8), intent(in) :: am(pcols,nrfm)               ! Mean air molecular weight interpolated to Fomichev grid
      real(r8), intent(in) :: co2(pcols,nrfm)              ! CO2 vmr interpolated to Fomichev grid
      real(r8), intent(in) :: djm(pcols,nrfm)              ! DJM coefficient in recurrence formula
      real(r8), intent(in) :: dj0(pcols,nrfm)              ! DJ0 coefficient in recurrence formula
      real(r8), intent(in) :: aajm(pcols,nrfm)             ! AAJM coefficient in recurrence formula
      real(r8), intent(in) :: aaj0(pcols,nrfm)             ! AAJ0 coefficient in recurrence formula
      real(r8), intent(in) :: alam(pcols,nrfm)             ! LAMBDA


!     Local variables
      real(r8) FU(pcols,nrfm)
      real(r8) FO3(pcols,nrfm)
      real(r8) H1(pcols)
      real(r8) H2(pcols)
      real(r8) H3(pcols)

      integer jj
      integer k
      integer i
      integer ks
      integer jjs


!     Output variables
      real(r8), intent(out) :: flux(pcols)                 ! Flux boundary condition for cool-to-space
      real(r8), intent(out) :: hco2(pcols,nrfm)            ! CO2 cooling in Fomichev grid
      real(r8), intent(out) :: ho3(pcols,nrfm)             ! O3 cooling in Fomichev grid

      real(r8) :: amat(nrfmlte,nrfmltelv)
      real(r8) :: bmat(nrfmlte,nrfmltelv)

!--------------------------------------------------------------------
!  update the amat and bmat matrices with time-dependent surf CO2 
!--------------------------------------------------------------------
      call set_matrices(amat,bmat)


      hco2(1:ncol,1:nrfm)=0.0_r8
      ho3(1:ncol,1:nrfm)=0.0_r8
      flux(1:ncol)=0.0_r8

!                                                                    
! grid levels for height integration   (p.s.h. distance = 0.25*IG)   
!

      do k=1,nrfm
         do i=1,ncol                                                 
            fu(i,k)=exp(-960.217_r8/tv(i,k))                                 
            fo3(i,k)=exp(-1500._r8/tv(i,k))                                  
         enddo
      enddo
     
!                                                                     
! calculate the heating rates for layer below s.h.p. = 12.5           
!   15 um CO2 + 9.6 um O3:                                            
!                                                                     
!                                                                     
!     **** COOLING RATE IN BOTH O3 AND CO2 BANDS (X=2-10.5) MATRIX    
!     **** APPROACH                                                   
!

!     Adding KS=K+8 maps NFRMC into NFRM
      do  k=1,5                                                      
         ks = k+8                                                     
         do i=1,ncol                                                  
            h2(i) = (amat(k,1)+bmat(k,1)*fu(i,ks))*fu(i,1)            
            h3(i) = ao3(k,1)*fo3(i,1)                                 
         enddo
         do jj=3,nrfmltelv
            jjs = ks+ig(jj)                                             
            do i=1,ncol                                              
               h2(i) = h2(i)+(amat(k,jj)+bmat(k,jj)*fu(i,ks))*fu(i,jjs)    
               h3(i) = h3(i)+ao3(k,jj)*fo3(i,jjs)                          
            enddo
         enddo
         do  i=1,ncol                                               
            hco2(i,ks) = h2(i)                                           
            ho3(i,ks) = h3(i)*o3(i,ks)                                   
         enddo
      enddo

      do k=6,18                                                      
         ks = k+8                                                     
         do i=1,ncol                                               
            h2(i) = (amat(k,1)+bmat(k,1)*fu(i,ks))*fu(i,1)               
            h3(i) = ao3(k,1)*fo3(i,1)                                    
         enddo
         do jj=2,nrfmltelv
            jjs = ks+ig(jj)                                              
            do i=1,ncol                                               
               h2(i) = h2(i)+(amat(k,jj)+bmat(k,jj)*fu(i,ks))*fu(i,jjs)     
               h3(i) = h3(i)+ao3(k,jj)*fo3(i,jjs)                           
            enddo
         enddo
         do  i=1,ncol                                            
            hco2(i,ks) = h2(i)                                            
            ho3(i,ks) = h3(i)*o3(i,ks)                                    
         enddo
      enddo

      do k=19,35                                                        
         ks = k+8                                                        
         do i=1,ncol                                                  
            h2(i) = 0._r8                                                      
            h3(i) = 0._r8                                                      
         enddo
         do  jj=1,nrfmltelv
            jjs = ks+ig(jj)                                                 
            do  i=1,ncol                                                 
               h2(i) = h2(i)+(amat(k,jj)+bmat(k,jj)*fu(i,ks))*fu(i,jjs)        
               h3(i) = h3(i)+ao3(k,jj)*fo3(i,jjs)                              
            enddo
         enddo
         do  i=1,ncol                                                  
            hco2(i,ks) = h2(i)                                               
            ho3(i,ks) = h3(i)*o3(i,ks)                                       
         enddo
      enddo

!                                                                         
!     **** COOLING RATE IN CO2 BANDS (X=10.75-12.5, MATRIX APPROACH)      
!                                                                         
      do k=36,43 
                                                  
         ks = k+8                                                        

         do  i=1,ncol                                                 
            h2(i) = 0._r8                                                      
         enddo

         do  jj=1,nrfmltelv
            jjs = ks+ig(jj)
            do i=1,ncol                                                 
               h2(i) = h2(i) + ( amat(k,jj) + bmat(k,jj)*fu(i,ks) ) * fu(i,jjs)
            enddo
         enddo

         do  i=1,ncol
            hco2(i,ks) = h2(i)
            ho3(i,ks) = 0._r8
         enddo

      enddo


!     Define boundary condition at XR=12.5 for recurrence formula
      do k=1,nrfm
         if (xr(k).eq.12.5_r8) then
            do i=1,ncol
               h1(i)=hco2(i,k)/(co2(i,k)*(1._r8-alam(i,k))*constb)
            enddo
         endif
      enddo


!     Do the rest of the XR domain 
!     (transition region 12.75 <= p.s.h. <= 16.5)
      do k=1,nrfm

         if (xr(k).ge.12.75_r8 .and. xr(k).le.16.5_r8) then

            do i=1,ncol

               h2(i)  = ( aajm(i,k)*h1(i) + djm(i,k)*fu(i,k-1)    &
                    - dj0(i,k)*fu(i,k) ) / aaj0(i,k)

               hco2(i,k) = h2(i)*co2(i,k)*(1._r8-alam(i,k))/am(i,k)*const
               ho3(i,k)  = 0._r8

               h1(i)     = h2(i)

            enddo               ! next longitude
            
            if (xr(k).eq.16.5_r8) then

!     Calculate FLUX at the top of the transition region (XR=16.5)
               do  i=1,ncol
                  flux(i) = h2(i) + fu(i,nrfm)
               enddo

               
            endif

         endif

      enddo                     ! next NLTE level


      return
  end subroutine viccooln

!============================================================================================

      subroutine cool2space (ncol,t,mwair,vn2,vo2,vo,vco2,ndenair,xnorm,flux,hc2s)

!
!     Adapted from Ray Roble's model by F. Sassi (Nov. 1999)
!     
!     Performs cool-to-space cooling calculations
!     This mod operates on the same vertical grid of the GCM
!
!-----------------------------------------------------------------
implicit none
!-----------------------------------------------------------------


!     Input variables
      integer, intent(in) :: ncol                          ! number of atmospheric columns

      real(r8), intent(in) :: t(pcols,pver)                       ! neutral temperature
      real(r8), intent(in) :: vn2(pcols,pver)                     ! N2 vmr
      real(r8), intent(in) :: vo2(pcols,pver)                     ! O2 vmr
      real(r8), intent(in) :: vco2(pcols,pver)                    ! CO2 vmr
      real(r8), intent(in) :: vo(pcols,pver)                      ! O vmr
      real(r8), intent(in) :: ndenair(pcols,pver)                 ! mean air no. density
      real(r8), intent(in) :: mwair(pcols,pver)                   ! mean air molecular weight
      real(r8), intent(in) :: flux(pcols)                         ! Radiative flux at the top of the NLTE region
      real(r8), intent(in) :: xnorm(pcols,pver)                        ! p.s.h.

!     Output  variables
      real(r8), intent(out) :: hc2s(pcols,pver)                   ! cool-to-space cooling

!     Local variables
      real(r8) tt
      real(r8) y
      real(r8) zn2
      real(r8) zo2
      real(r8) zz
      real(r8) alam
      real(r8) rko

      integer k
      integer i      

!     ****                                        
!     ****   CO2 COOL-TO-SPACE APPROXIMATION      
!     **** 


      hc2s(1:ncol,1:pver)=0.0_r8
      do k = 1,pver

         do i=1,ncol

            if (xnorm(i,k) .gt. 16.5_r8) then

               tt     = t(i,k)                                        
               y      = tt**(-1._r8/3._r8)                                  
               zn2    = 5.5e-17_r8*sqrt(tt) + 6.7e-10_r8 * exp(-83.8_r8*y)    
               zo2    = 1.0e-15_r8*exp(23.37_r8 - 230.9_r8*y + 564._r8*y*y)
               rko=3.0e-12_r8 

!     
!     ****  COLLISIONAL DEACTIVATION RATE:                       
!
               zz     = (vn2(i,k)*zn2 + vo2(i,k)*zo2 + vo (i,k)*rko)*ndenair(i,k)

!                                                                
!     ****                                                       
!                                                                
               alam   = a10/(a10+zz)
               hc2s(i,k) = const/mwair(i,k)*vco2(i,k)*(1._r8-alam)      &
               *(flux(i)-exp(-960.217_r8/tt))

            endif
         
         enddo                  ! end-loop in longitude

      enddo                     ! end-loop in levels

      return
 end subroutine cool2space


!=====================================================================


real(r8) function a18lin (x,xn,yn,m,n)

! input:
!  X - argument for which a value of function should be found
!  XN(N),YN(N) - values of function YN(N) at XN(N) grid. X(N) should be
!                ordered so that X(I-1) < X(I).
! output:
!  A18LIN - value of function for X

      implicit none

!
! Args:
      integer,intent(in) :: m,n
      real(r8),intent(in) :: x
      real(r8),intent(in) :: xn(n)
      real(r8),intent(in) :: yn(n)
!
! Local:
      integer :: k,i

      k=m-1
      the_loop: do i=m,n
         k=k+1
         if (x-xn(i).le.0._r8) exit the_loop
      enddo the_loop
      if(k.eq.1) k=2

! k has been found so that xn(k).le.x.lt.xn(k+1)

      a18lin=(yn(k)-yn(k-1))/(xn(k)-xn(k-1))*(x-xn(k))+yn(k)

     return

   end function a18lin

!=============================================================================


      subroutine a18int(x1,y1,x2,y2,n1,n2)

!
! third order spline interpolation
! input argument and function:  X1(1:N1),Y1(1:N1)
! output argument and function: X2(1:N2)X2(1:N2),Y2(1:N2)
! the necessary conditionts are: X1(I) < X1(I+1), and the same for X2 array.
!

      implicit none
!
! Args:
      integer,  intent(in) :: n1
      integer,  intent(in) :: n2
      real(r8), intent(in) :: x1(n1)
      real(r8), intent(in) :: y1(n1)
      real(r8), intent(in) :: x2
      real(r8), intent(out) :: y2
!
! Local:
      real(r8) :: a(150),e(150),f(150),h(150),h2,h1,f1,f2,f3,g
      integer :: nvs,k,kr,l
!
      h2=x1(1)
      nvs=n1-1
      do 1 k=1,nvs
      h1=h2
      h2=x1(k+1)
      h(k)=h2-h1
    1 continue
      a(1)=0._r8
      a(n1)=0._r8
      e(n1)=0._r8
      f(n1)=0._r8
      h1=h(n1-1)
      f1=y1(n1-1)
      f2=y1(n1)
      do 2 kr=2,nvs
      k=nvs+2-kr
      h2=h1
      h1=h(k-1)
      f3=f2
      f2=f1
      f1=y1(k-1)
      g=1._r8/(h2*e(k+1)+2._r8*(h1+h2))
      e(k)=-h1*g
      f(k)=(3._r8*((f3-f2)/h2-(f2-f1)/h1)-h2*f(k+1))*g
    2 continue
      g=0._r8
      do 3 k=2,nvs
      g=e(k)*g+f(k)
      a(k)=g
    3 continue
      l=1
!$$$      do 4 i=1,n2
!$$$      g=x2(i)
      g=x2
      do 6 k=l,nvs
      if(g.gt.x1(k+1)) goto 6
      l=k
      goto 5
    6 continue
      l=nvs
    5 g=g-x1(l)
      h2=h(l)
      f2=y1(l)
      f1=h2**2
      f3=g**2
      y2=f2+g/h2*(y1(l+1)-f2-(a(l+1)*(f1-f3)+         &
                a(l)*(2._r8*f1-3._r8*g*h2+f3))/3._r8)
!$$$    4 continue
      return
   end  subroutine a18int

!==================================================================================================

   subroutine nocooling(ncol,t,pmid,nommr,o1mmr,o2mmr,o3mmr,n2mmr,nocool)

!
! Calculate NO cooling (ref: Kockarts, GRL, vol. 7, pp 137-140, 1980)
!

     integer, intent(in) :: ncol                  ! number of column in chunck

     real(r8), intent(in) :: t(pcols,pver)        ! neutral gas temperature (K)
     real(r8), intent(in) :: pmid(pcols,pver)     ! model pressure at midpoints (Pa)
     real(r8), intent(in) :: nommr(pcols,pver)    ! NO (in mmr)
     real(r8), intent(in) :: o1mmr(pcols,pver)    ! O (in mmr)
     real(r8), intent(in) :: o2mmr(pcols,pver)    ! O2 (in mmr)
     real(r8), intent(in) :: o3mmr(pcols,pver)    ! O3 (in mmr)
     real(r8), intent(in) :: n2mmr(pcols,pver)    ! N2 (in mmr)

     real(r8), intent(out) :: nocool(pcols,pver)  ! NO-cooling (K/S)

! Local space
     real(r8) :: mwair(pcols,pver)                ! mean molecular weight
     real(r8) :: pres_cgs(pcols,pver)             ! pressure in cgs units
     real(r8) :: nair(pcols,pver)                 ! mean air number density (molecules/cm3)
     real(r8) :: o1vmr(pcols,pver)                ! O (vmr)
     real(r8) :: o2vmr(pcols,pver)                ! O2 (vmr)
     real(r8) :: no_conc                          ! NO concentration divided by mean air density
     real(r8) :: no_deact                         ! effective NO deactivation rate multiplied by air concentration

     ! Hwang et al., "Vibrational relaxation of NO(v=1) by oxygen atoms between 295 and 825 K"
     ! JGR, Vol. 108, No. A3, 1109, doi:10.1029/2002JA009688, 2003
     real(r8), parameter :: o1_rate = 4.2e-11_r8     ! O1 reaction rate (cm3/sec)
     real(r8), parameter :: o2_rate = 2.4e-14_r8     ! O2 reaction rate (cm3/sec)
     real(r8), parameter :: phot_e  = 3.726e-13_r8   ! photon energy at 5.3 mum (erg)
     real(r8), parameter :: trans_prob = 13.3_r8     ! Einstein transition probability

     integer :: i,k
!---------------------------------------------------------------------------------

! Calculate mean air molecular weight
     do k=1, pver
        do i=1, ncol
           mwair(i,k) = 1._r8 /                                                                 &
                (o2mmr(i,k) /o2_mw + (o1mmr(i,k)+o3mmr(i,k))/o1_mw + n2mmr(i,k)/n2_mw)
        end do
     end do

     do k=1,pver
        do i=1,ncol
! convert mmr to vmr
           o1vmr(i,k) = o1mmr(i,k) * mwair(i,k) / o1_mw
           o2vmr(i,k) = o2mmr(i,k) * mwair(i,k) / o2_mw
! Convert pressure to cgs units
           pres_cgs(i,k) = pmid(i,k) * 10._r8
! calculate mean air number density
           nair(i,k) = anav * pres_cgs(i,k) / (ur*t(i,k)) 
        end do
     end do

!----------------------------------------------------------------------------------
! NO cooling
!----------------------------------------------------------------------------------
     do k=1,pver
        do i=1,ncol
!----------------------------------------------------------------------------------
! calcualte effective NO deactivation rate times mean air concentration
!----------------------------------------------------------------------------------
           no_deact =  nair(i,k) * (o1_rate * o1vmr(i,k) + o2_rate * o2vmr(i,k) )
!----------------------------------------------------------------------------------
! calculate NO concentration:
!
!                                   Na
!                n(NO) = mmr(NO) * ------ * rho(air)
!                                   M(NO)
!
! where M(NO) is NO molecular weight and rho(air) is mean air density.
! However, in order to produce a heating in energy per unit mass, we need
! to divide by rho(air) at the end. Therefore, rho(air) IS NOT INCLUDED IN
! THE CALCULATION OF n(NO).
!----------------------------------------------------------------------------------
           no_conc = anav * nommr(i,k) / no_mw 
!----------------------------------------------------------------------------------
! Heating calculation. It produces a heating in erg/g/s.
! convert to J/kg/s
!----------------------------------------------------------------------------------
           nocool(i,k) = -1.e-4_r8 * phot_e * trans_prob * no_conc &
                                  * (no_deact / (no_deact + trans_prob)) * exp(-2700._r8/t(i,k))
        enddo
     enddo
     

!
   end subroutine nocooling

!==================================================================================================

   subroutine o3pcooling( ncol, t, o1mmr, o3pcool )
     !
     ! Adapted from TIE-GCM 
     ! Original equation is from Bates [1951]
     ! Masking factors are from Kockarts and Peetermans [1981]
     !
     
     integer, intent(in) :: ncol                  ! number of column in chunck

     real(r8), intent(in) :: t(pcols,pver)        ! neutral gas temperature (K)
     real(r8), intent(in) :: o1mmr(pcols,pver)    ! O (in mmr)

     real(r8), intent(out) :: o3pcool(pcols,pver)  ! O(3p)-cooling (K/S)

     ! Local space

     integer  :: i,k
     real(r8) :: invtemp(ncol,pver)
     real(r8) :: work1(ncol,pver)             
     real(r8) :: work2(ncol,pver)            
     real(r8) :: anavfac

     real(r8),parameter :: &
          an(3) = (/0.835E-18_r8, 0.6_r8, 0.2_r8/), &  
          bn(3) = (/228._r8,228._r8,325._r8/)      ! coefficients in Bates equation

     invtemp(:ncol,:) = 1.0_r8/t(:ncol,:)
     anavfac = anav/o1_mw

     do k=1,pver
        do i=1,ncol
           work1(i,k) = an(1)*o3pxfac(k)*anavfac
           work1(i,k) = work1(i,k)*o1mmr(i,k)*exp(-bn(1)*invtemp(i,k))
           work2(i,k) = 1._r8 + an(2)*exp(-bn(2)*invtemp(i,k)) &
                              + an(3)*exp(-bn(3)*invtemp(i,k))
           o3pcool(i,k) = -work1(i,k)/work2(i,k)   ! erg/g/s
           o3pcool(i,k) = o3pcool(i,k) * 1E-04_r8  ! convert units from erg/g/s to J/kg/s
        enddo
     enddo

   end subroutine o3pcooling

end module nlte_fomichev
